In 2017, the Tuveson Lab at Cold Spring Harbor Cancer Center published a paper written by Elyada et al (2017) that detailed the discovery of cancer-associated fibroblasts (CAFs) in mice. The subtypes were then validated in human samples affected with PDAC in a subsequent paper released in 2019. Here we will use SCISSORS to identify the CAF subtypes within the larger fibroblast population, immune cell types within broadly-defined immune clusters, and ductal & PDAC subtypes within the ductal group.
library(pals) # high N discrete color palettes
library(VAM) # single cell GSEA
library(dplyr) # tidy data
library(Seurat) # single cell infrastructure
library(ggplot2) # plots
library(SingleR) # cell type assignment
library(decoderr) # de novo deconvolution
library(SCISSORS) # our package
library(mixtools) # Gaussian mixture model estimation
library(patchwork) # plot grids
library(paletteer) # more color palettes
library(CONICSmat) # CNV estimation
library(reticulate) # Python interface
library(wesanderson) # even more color palettesimport numpy as np
from openTSNE import TSNEEmbedding
from openTSNE import initialization
from openTSNE.affinity import Multiscale
from openTSNE.affinity import PerplexityBasedNNFirst we load in the \(\text{gene} \times \text{cell}\) counts matrix, then create a Seurat object to hold it in. Next, we add samplename, tissue type, and patient sex metadata taken from the publicly available dataset.
raw_counts <- Read10X(data.dir = "~/Desktop/Data/Elyada Raw/All Human/")
pdac <- CreateSeuratObject(raw_counts,
project = "Elyada",
min.cells = 3,
min.features = 500)
pdac@meta.data$sample <- case_when(grepl("-1", rownames(pdac@meta.data)) ~ "SRR9274536",
grepl("-2", rownames(pdac@meta.data)) ~ "SRR9274537",
grepl("-3", rownames(pdac@meta.data)) ~ "SRR9274538",
grepl("-4", rownames(pdac@meta.data)) ~ "SRR9274539",
grepl("-5", rownames(pdac@meta.data)) ~ "SRR9274540",
grepl("-6", rownames(pdac@meta.data)) ~ "SRR9274541",
grepl("-7", rownames(pdac@meta.data)) ~ "SRR9274542",
grepl("-8", rownames(pdac@meta.data)) ~ "SRR9274543",
grepl("-9", rownames(pdac@meta.data)) ~ "SRR9274544")
pdac@meta.data$condition <- case_when(grepl("-1", rownames(pdac@meta.data)) ~ "PDAC",
grepl("-2", rownames(pdac@meta.data)) ~ "PDAC",
grepl("-3", rownames(pdac@meta.data)) ~ "PDAC",
grepl("-4", rownames(pdac@meta.data)) ~ "PDAC",
grepl("-5", rownames(pdac@meta.data)) ~ "PDAC",
grepl("-6", rownames(pdac@meta.data)) ~ "AdjNorm",
grepl("-7", rownames(pdac@meta.data)) ~ "PDAC",
grepl("-8", rownames(pdac@meta.data)) ~ "AdjNorm",
grepl("-9", rownames(pdac@meta.data)) ~ "PDAC")
pdac@meta.data$sex <- case_when(grepl("-1", rownames(pdac@meta.data)) ~ "female",
grepl("-2", rownames(pdac@meta.data)) ~ "male",
grepl("-3", rownames(pdac@meta.data)) ~ "male",
grepl("-4", rownames(pdac@meta.data)) ~ "male",
grepl("-5", rownames(pdac@meta.data)) ~ "male",
grepl("-6", rownames(pdac@meta.data)) ~ "female",
grepl("-7", rownames(pdac@meta.data)) ~ "female",
grepl("-8", rownames(pdac@meta.data)) ~ "male",
grepl("-9", rownames(pdac@meta.data)) ~ "female")We run the typical single cell pre-processing steps on our cells - normalization, dimension reduction, and clustering. The t-SNE embedding looks decent, but could definitely be improved. Globally, the clusters are arranged in a blob - which isn’t very informative - though the local structure seems to have been preserved fairly well. Clusters 5 and 11 are split in two, which we would prefer not to have happen.
pdac <- PrepareData(pdac,
n.variable.genes = 4000,
n.PC = 20,
which.dim.reduc = "tsne",
initial.resolution = .5,
random.seed = 629)
p0 <- DimPlot(pdac, reduction = "tsne") +
scale_color_manual(values = unname(tableau20())) +
labs(x = "t-SNE 1", y = "t-SNE 2") +
theme_yehlab() +
guides(color = guide_legend(nrow = 2, override.aes = list(size = 3)))
p0## [1] "Normalizing counts using SCTransform"
## [1] "Running t-SNE on 20 principal components with perplexity = 30"
## [1] "Clustering cells in PCA space using k ~ 155 & resolution = 0.5"
## [1] "Found 13 unique clusters"
I think the two-dimensional visualization of the cells could be improved. We’ll try using UMAP and the Fast Fourier Transform-accelerated Fit-SNE (as implemented in the openTSNE library) to improve the embedding.
It seems like UMAP does a good job of clearly separating our clusters and preserving the global structure of the data. However, difficult to see local structure within some of the clusters due to their density, and cluster 5 is split across the 2D UMAP space.
pdac <- RunUMAP(pdac,
reduction = "pca",
dims = 1:20,
umap.method = "uwot",
n.components = 2,
n.epochs = 750,
n.neighbors = 50,
metric = "cosine",
seed.use = 629,
verbose = FALSE)
p1 <- DimPlot(pdac, reduction = "umap") +
scale_color_manual(values = unname(tableau20())) +
labs(x = "UMAP 1", y = "UMAP 2") +
theme_yehlab() +
guides(color = guide_legend(nrow = 2, override.aes = list(size = 3)))
p1For information on how to install the openTSNE implementation of Fit-SNE and how to run the algorithm, please visit the excellent GitHub repository of Pavlin Policar.
First we’ll run a simple, standard Fit-SNE embedding. It’s necessary to make the PCA embeddings accessible by Python.
pc_df <- Embeddings(pdac, reduction = "pca")# import data
pc_df = np.array(r.pc_df)
# run Fit-SNE
affin = PerplexityBasedNN(pc_df, perplexity=50, metric='cosine', random_state=629)
init = initialization.pca(pc_df, random_state=629)
tsne1 = TSNEEmbedding(init, affin, negative_gradient_method='fft')
embed1 = tsne1.optimize(n_iter=350, exaggeration=12, momentum=0.6)
embed2 = embed1.optimize(n_iter=1000, momentum=0.8)The embedding looks good, but the subclusters within clusters 5 and 11 are still clearly separated, and a small part of cluster 12 is far from the rest of the cells in the cluster.
embed <- as.matrix(py$embed2)
rownames(embed) <- colnames(pdac)
pdac@reductions$fitsne <- CreateDimReducObject(embeddings = embed,
key = "FitSNE_",
assay = "SCT",
global = TRUE)
p2 <- DimPlot(pdac, reduction = "fitsne") +
scale_color_manual(values = unname(tableau20())) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
guides(color = guide_legend(nrow = 2, override.aes = list(size = 3)))
p2We use the ability of Fit-SNE to consider multiple perplexities in an attempt to provide a more accurate representation of both local and global structure in our data.
affin_anneal = PerplexityBasedNN(pc_df, perplexity=150, metric='cosine', random_state=629)
tsne2 = TSNEEmbedding(init, affin_anneal, negative_gradient_method='fft')
embed3 = tsne2.optimize(n_iter=350, exaggeration=12, momentum=0.5)
embed4 = embed1.optimize(n_iter=1000, exaggeration=1, momentum=0.8)
affin_anneal.set_perplexity(30)
embed5 = embed4.optimize(n_iter=500, momentum=0.8)The perplexity annealing hasn’t really improved or even changed the embedding much at all.
embed_multi <- as.matrix(py$embed5)
rownames(embed_multi) <- colnames(pdac)
pdac@reductions$fitsne_multi <- CreateDimReducObject(embeddings = embed_multi,
key = "FitSNE_",
assay = "SCT",
global = TRUE)
p3 <- DimPlot(pdac, reduction = "fitsne_multi") +
scale_color_manual(values = unname(tableau20())) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
guides(color = guide_legend(nrow = 2, override.aes = list(size = 3)))
p3Finally, we’ll use a multiscale kernel to represent the data in the high-dimensional space instead of the default Gaussian. This in tandem with the multiple perplexity values will preserve as much as possible of both local and global structure.
affin_multi = Multiscale(pc_df, perplexities=[30, 200], metric='cosine', random_state=629)
tsne3 = TSNEEmbedding(init, affin_multi, negative_gradient_method='fft')
embed6 = tsne3.optimize(n_iter=450, exaggeration=12, momentum=0.6)
embed7 = embed6.optimize(n_iter=1200, exaggeration=1, momentum=0.8)The multiscale trick doesn’t completely fix the embedding issues either, but it does clean up a bit of the noise in the plot, and it does place the subclusters in clusters 11 and 5 closer together. We’ll use this embedding going forward and make sure to examine those clusters for subgroups using SCISSORS.
embed_kern <- as.matrix(py$embed7)
rownames(embed_kern) <- colnames(pdac)
pdac@reductions$fitsne_kern<- CreateDimReducObject(embeddings = embed_kern,
key = "FitSNE_",
assay = "SCT",
global = TRUE)
p4 <- DimPlot(pdac, reduction = "fitsne_kern") +
scale_color_manual(values = unname(tableau20())) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
guides(color = guide_legend(nrow = 1, override.aes = list(size = 3)))
p4Looking at the embedding, we can see that clusters 3, 4, 5, 10, and 11 contain clear small subclusters. I’m also intrigued by the shape of cluster 8. We’ll be sure to investigate each of these in turn, and annotate the groups using marker genes.
Due to the way Seurat accesses cell embeddings, we’ll need to replace our original t-SNE dimension reduction in our Seurat object with the new Fit-SNE version. We’ll keep the original Barnes-Hut t-SNE embedding under a separate name.
pdac@reductions$bh_tsne <- pdac@reductions$tsne
pdac@reductions$tsne <- pdac@reductions$fitsne_kernHere we use the SingleR package to identify broad cell types. The reference dataset we load is an sctransform-normalized version of the raw counts available in scRNAseq::BaronPancreasData(), which consists of normal pancreas cells that were sequenced and annotated by the researchers.
sc_ref <- readRDS("/Volumes/labs/Home/Jen Jen Yeh Lab/Jack/scRNAseq/Seurat/single_cell_ref_normalized.Rds")
sc_preds <- SingleR(test = data.frame(pdac@assays$SCT@data),
ref = sc_ref,
labels = sc_ref$label,
method = "cluster",
clusters = pdac$seurat_clusters,
de.method = "wilcox")
pdac[["SingleR.labels.sc"]] <- sc_preds$labels[match(pdac[[]][["seurat_clusters"]], rownames(sc_preds))]
pdac$SingleR.labels.sc <- case_when(pdac$SingleR.labels.sc == "acinar" ~ "Acinar",
pdac$SingleR.labels.sc == "activated_stellate" ~ "Activated Stellate",
pdac$SingleR.labels.sc == "ductal" ~ "Ductal",
pdac$SingleR.labels.sc == "macrophage" ~ "Macrophage",
pdac$SingleR.labels.sc == "t_cell" ~ "T",
pdac$SingleR.labels.sc == "quiescent_stellate" ~ "Quiescent Stellate")We can see that there’s a large immune population, as well as smaller ductal, fibroblast (denoted activated and quiescent stellate in the reference dataset), and acinar groups. These broad cell types line up with what we expected to see given the authors’ original cell cluster annotations.
p5 <- DimPlot(pdac, reduction = "tsne", group.by = "SingleR.labels.sc") +
scale_color_manual(values = paletteer_d("miscpalettes::brightPastel")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme(plot.title = element_blank()) +
theme_yehlab() +
guides(color = guide_legend(nrow = 2, override.aes = list(size = 3)))
p5This dataset is composed of labeled & log-normalized bulk RNA-seq samples from the Human Primary Cell Atlas.
bulk_ref <- HumanPrimaryCellAtlasData()
bulk_preds <- SingleR(test = data.frame(pdac@assays$SCT@data),
ref = bulk_ref,
labels = bulk_ref$label.main,
method = "cluster",
clusters = pdac$seurat_clusters,
de.method = "wilcox")
pdac[["SingleR.labels.bulk"]] <- bulk_preds$labels[match(pdac[[]][["seurat_clusters"]],
rownames(bulk_preds))]
pdac$SingleR.labels.bulk <- case_when(pdac$SingleR.labels.bulk == "B_cell" ~ "B",
pdac$SingleR.labels.bulk == "Epithelial_cells" ~ "Epithelial",
pdac$SingleR.labels.bulk == "B_cell-" ~ "B",
pdac$SingleR.labels.bulk == "T_cells" ~ "T",
pdac$SingleR.labels.bulk == "Endothelial_cells" ~ "Endothelial",
TRUE ~ pdac$SingleR.labels.bulk)The bulk reference gives us somewhat more granular labels for the immune cells, and confirms the identities of the ductal / epithelial and stroma clusters. We see that the clusters denoted quiescent stellate by the single cell reference are here labeled as endothelial cells, and the activated stellate cells are confirmed to be fibroblasts.
p6 <- DimPlot(pdac, reduction = "tsne", group.by = "SingleR.labels.bulk") +
scale_color_manual(values = paletteer_d("miscpalettes::brightPastel")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme(plot.title = element_blank()) +
theme_yehlab() +
guides(color = guide_legend(nrow = 2, override.aes = list(size = 3)))
p6One shouldn’t use SingleR as the final authority for cell types, but we were able to confirm the identities of the ductal and fibroblast clusters, which is important for this dataset as the cells we are most interested in are the cancer-associated fibroblasts (CAFs).
Next we’ll attempt to identify malignant cells using single-cell copy number variation estimation as implemented in the CONCISmat package. Details of the GMM methodology used can be found at the Diaz Lab’s GitHub repository. Note: this step is memory-intensive because 1) it requires the sparse counts matrix to be cast to a dense matrix and 2) a lot of Gaussian mixture models get estimated. If your machine doesn’t have a lot of RAM it might be best to skip this and manually annotate the malignant cells instead through the usage of known canonical markers or the VAM single-cell GSEA methodology.
chrom_regions <- read.table("/Volumes/labs/Home/Jen Jen Yeh Lab/Jack/scRNAseq/chrom_arm_positions.txt",
sep = "\t",
row.names = 1,
header = TRUE)
gene_pos <- getGenePositions(rownames(pdac))
cpm <- t(t(as.matrix(pdac@assays$SCT@counts)) / colSums(as.matrix(pdac@assays$SCT@counts))) * 10^5
cpm <- log2(cpm + 1)
norm_factor <- calcNormFactors(cpm)
cnv_est <- plotAll(mat = cpm,
normFactor = norm_factor,
regions = chrom_regions,
gene_pos = gene_pos,
fname = "Elyada")## [1] "Fitting GMM for chr1 0:122026459 iteration 1"
## number of iterations= 404
## [1] "Fitting GMM for chr1 0:122026459 iteration 2"
## number of iterations= 561
## [1] "Fitting GMM for chr1 0:122026459 iteration 3"
## number of iterations= 437
## [1] "Fitting GMM for chr1 0:122026459 iteration 4"
## number of iterations= 493
## [1] 44
## [1] "Fitting GMM for chr1 124932724:248956422 iteration 1"
## number of iterations= 97
## [1] "Fitting GMM for chr1 124932724:248956422 iteration 2"
## number of iterations= 581
## [1] "Fitting GMM for chr1 124932724:248956422 iteration 3"
## number of iterations= 240
## [1] "Fitting GMM for chr1 124932724:248956422 iteration 4"
## number of iterations= 515
## [1] 44
## [1] "Fitting GMM for chr2 0:92188145 iteration 1"
## WARNING! NOT CONVERGENT!
## number of iterations= 1000
## [1] "Fitting GMM for chr2 0:92188145 iteration 2"
## WARNING! NOT CONVERGENT!
## number of iterations= 1000
## [1] "Fitting GMM for chr2 0:92188145 iteration 3"
## WARNING! NOT CONVERGENT!
## number of iterations= 1000
## [1] "Fitting GMM for chr2 0:92188145 iteration 4"
## WARNING! NOT CONVERGENT!
## number of iterations= 1000
## [1] 44
## [1] "Fitting GMM for chr2 94090557:242193529 iteration 1"
## number of iterations= 786
## [1] "Fitting GMM for chr2 94090557:242193529 iteration 2"
## number of iterations= 869
## [1] "Fitting GMM for chr2 94090557:242193529 iteration 3"
## WARNING! NOT CONVERGENT!
## number of iterations= 1000
## [1] "Fitting GMM for chr2 94090557:242193529 iteration 4"
## WARNING! NOT CONVERGENT!
## number of iterations= 1000
## [1] 44
## [1] "Fitting GMM for chr3 0:90772458 iteration 1"
## number of iterations= 651
## [1] "Fitting GMM for chr3 0:90772458 iteration 2"
## number of iterations= 763
## [1] "Fitting GMM for chr3 0:90772458 iteration 3"
## number of iterations= 697
## [1] "Fitting GMM for chr3 0:90772458 iteration 4"
## number of iterations= 674
## [1] 44
## [1] "Fitting GMM for chr3 93655574:198295559 iteration 1"
## number of iterations= 149
## [1] "Fitting GMM for chr3 93655574:198295559 iteration 2"
## number of iterations= 149
## [1] "Fitting GMM for chr3 93655574:198295559 iteration 3"
## number of iterations= 174
## [1] "Fitting GMM for chr3 93655574:198295559 iteration 4"
## number of iterations= 144
## [1] 44
## [1] "Fitting GMM for chr4 0:49712061 iteration 1"
## WARNING! NOT CONVERGENT!
## number of iterations= 1000
## [1] "Fitting GMM for chr4 0:49712061 iteration 2"
## WARNING! NOT CONVERGENT!
## number of iterations= 1000
## [1] "Fitting GMM for chr4 0:49712061 iteration 3"
## WARNING! NOT CONVERGENT!
## number of iterations= 1000
## [1] "Fitting GMM for chr4 0:49712061 iteration 4"
## number of iterations= 880
## [1] 44
## [1] "Fitting GMM for chr4 51743951:190214555 iteration 1"
## number of iterations= 427
## [1] "Fitting GMM for chr4 51743951:190214555 iteration 2"
## number of iterations= 316
## [1] "Fitting GMM for chr4 51743951:190214555 iteration 3"
## number of iterations= 336
## [1] "Fitting GMM for chr4 51743951:190214555 iteration 4"
## number of iterations= 415
## [1] 44
## [1] "Fitting GMM for chr5 0:46485900 iteration 1"
## WARNING! NOT CONVERGENT!
## number of iterations= 1000
## [1] "Fitting GMM for chr5 0:46485900 iteration 2"
## WARNING! NOT CONVERGENT!
## number of iterations= 1000
## [1] "Fitting GMM for chr5 0:46485900 iteration 3"
## WARNING! NOT CONVERGENT!
## number of iterations= 1000
## [1] "Fitting GMM for chr5 0:46485900 iteration 4"
## number of iterations= 935
## [1] 44
## [1] "Fitting GMM for chr5 50059807:181538259 iteration 1"
## number of iterations= 411
## [1] "Fitting GMM for chr5 50059807:181538259 iteration 2"
## number of iterations= 597
## [1] "Fitting GMM for chr5 50059807:181538259 iteration 3"
## number of iterations= 533
## [1] "Fitting GMM for chr5 50059807:181538259 iteration 4"
## number of iterations= 509
## [1] 44
## [1] "Fitting GMM for chr6 0:58553888 iteration 1"
## number of iterations= 404
## [1] "Fitting GMM for chr6 0:58553888 iteration 2"
## number of iterations= 363
## [1] "Fitting GMM for chr6 0:58553888 iteration 3"
## number of iterations= 788
## [1] "Fitting GMM for chr6 0:58553888 iteration 4"
## number of iterations= 563
## [1] 44
## [1] "Fitting GMM for chr6 59829934:170805979 iteration 1"
## WARNING! NOT CONVERGENT!
## number of iterations= 1000
## [1] "Fitting GMM for chr6 59829934:170805979 iteration 2"
## WARNING! NOT CONVERGENT!
## number of iterations= 1000
## [1] "Fitting GMM for chr6 59829934:170805979 iteration 3"
## WARNING! NOT CONVERGENT!
## number of iterations= 1000
## [1] "Fitting GMM for chr6 59829934:170805979 iteration 4"
## number of iterations= 767
## [1] 44
## [1] "Fitting GMM for chr7 0:58169653 iteration 1"
## WARNING! NOT CONVERGENT!
## number of iterations= 1000
## [1] "Fitting GMM for chr7 0:58169653 iteration 2"
## WARNING! NOT CONVERGENT!
## number of iterations= 1000
## [1] "Fitting GMM for chr7 0:58169653 iteration 3"
## WARNING! NOT CONVERGENT!
## number of iterations= 1000
## [1] "Fitting GMM for chr7 0:58169653 iteration 4"
## WARNING! NOT CONVERGENT!
## number of iterations= 1000
## [1] 44
## [1] "Fitting GMM for chr7 61528020:159345973 iteration 1"
## number of iterations= 580
## [1] "Fitting GMM for chr7 61528020:159345973 iteration 2"
## number of iterations= 507
## [1] "Fitting GMM for chr7 61528020:159345973 iteration 3"
## number of iterations= 579
## [1] "Fitting GMM for chr7 61528020:159345973 iteration 4"
## number of iterations= 567
## [1] 44
## [1] "Fitting GMM for chr8 0:44033744 iteration 1"
## number of iterations= 681
## [1] "Fitting GMM for chr8 0:44033744 iteration 2"
## number of iterations= 845
## [1] "Fitting GMM for chr8 0:44033744 iteration 3"
## number of iterations= 782
## [1] "Fitting GMM for chr8 0:44033744 iteration 4"
## number of iterations= 787
## [1] 44
## [1] "Fitting GMM for chr8 45877265:145138636 iteration 1"
## number of iterations= 718
## [1] "Fitting GMM for chr8 45877265:145138636 iteration 2"
## number of iterations= 614
## [1] "Fitting GMM for chr8 45877265:145138636 iteration 3"
## number of iterations= 614
## [1] "Fitting GMM for chr8 45877265:145138636 iteration 4"
## number of iterations= 614
## [1] 44
## [1] "Fitting GMM for chr9 0:43389635 iteration 1"
## WARNING! NOT CONVERGENT!
## number of iterations= 1000
## [1] "Fitting GMM for chr9 0:43389635 iteration 2"
## WARNING! NOT CONVERGENT!
## number of iterations= 1000
## [1] "Fitting GMM for chr9 0:43389635 iteration 3"
## WARNING! NOT CONVERGENT!
## number of iterations= 1000
## [1] "Fitting GMM for chr9 0:43389635 iteration 4"
## number of iterations= 863
## [1] 44
## [1] "Fitting GMM for chr9 45518558:138394717 iteration 1"
## number of iterations= 655
## [1] "Fitting GMM for chr9 45518558:138394717 iteration 2"
## number of iterations= 656
## [1] "Fitting GMM for chr9 45518558:138394717 iteration 3"
## number of iterations= 642
## [1] "Fitting GMM for chr9 45518558:138394717 iteration 4"
## number of iterations= 476
## [1] 44
## [1] "Fitting GMM for chr10 0:39686682 iteration 1"
## WARNING! NOT CONVERGENT!
## number of iterations= 1000
## [1] "Fitting GMM for chr10 0:39686682 iteration 2"
## WARNING! NOT CONVERGENT!
## number of iterations= 1000
## [1] "Fitting GMM for chr10 0:39686682 iteration 3"
## WARNING! NOT CONVERGENT!
## number of iterations= 1000
## [1] "Fitting GMM for chr10 0:39686682 iteration 4"
## WARNING! NOT CONVERGENT!
## number of iterations= 1000
## [1] 44
## [1] "Fitting GMM for chr10 41593521:133797422 iteration 1"
## number of iterations= 588
## [1] "Fitting GMM for chr10 41593521:133797422 iteration 2"
## number of iterations= 525
## [1] "Fitting GMM for chr10 41593521:133797422 iteration 3"
## number of iterations= 640
## [1] "Fitting GMM for chr10 41593521:133797422 iteration 4"
## number of iterations= 555
## [1] 44
## [1] "Fitting GMM for chr11 0:51078348 iteration 1"
## number of iterations= 236
## [1] "Fitting GMM for chr11 0:51078348 iteration 2"
## number of iterations= 293
## [1] "Fitting GMM for chr11 0:51078348 iteration 3"
## number of iterations= 294
## [1] "Fitting GMM for chr11 0:51078348 iteration 4"
## number of iterations= 158
## [1] 44
## [1] "Fitting GMM for chr11 54425074:135086622 iteration 1"
## number of iterations= 977
## [1] "Fitting GMM for chr11 54425074:135086622 iteration 2"
## number of iterations= 862
## [1] "Fitting GMM for chr11 54425074:135086622 iteration 3"
## number of iterations= 870
## [1] "Fitting GMM for chr11 54425074:135086622 iteration 4"
## number of iterations= 892
## [1] 44
## [1] "Fitting GMM for chr12 0:34769407 iteration 1"
## WARNING! NOT CONVERGENT!
## number of iterations= 1000
## [1] "Fitting GMM for chr12 0:34769407 iteration 2"
## WARNING! NOT CONVERGENT!
## number of iterations= 1000
## [1] "Fitting GMM for chr12 0:34769407 iteration 3"
## WARNING! NOT CONVERGENT!
## number of iterations= 1000
## [1] "Fitting GMM for chr12 0:34769407 iteration 4"
## WARNING! NOT CONVERGENT!
## number of iterations= 1000
## [1] 44
## [1] "Fitting GMM for chr12 37185252:133275309 iteration 1"
## number of iterations= 887
## [1] "Fitting GMM for chr12 37185252:133275309 iteration 2"
## number of iterations= 768
## [1] "Fitting GMM for chr12 37185252:133275309 iteration 3"
## number of iterations= 969
## [1] "Fitting GMM for chr12 37185252:133275309 iteration 4"
## number of iterations= 878
## [1] 44
## [1] "Fitting GMM for chr13 18051248:114364328 iteration 1"
## WARNING! NOT CONVERGENT!
## number of iterations= 1000
## [1] "Fitting GMM for chr13 18051248:114364328 iteration 2"
## WARNING! NOT CONVERGENT!
## number of iterations= 1000
## [1] "Fitting GMM for chr13 18051248:114364328 iteration 3"
## WARNING! NOT CONVERGENT!
## number of iterations= 1000
## [1] "Fitting GMM for chr13 18051248:114364328 iteration 4"
## WARNING! NOT CONVERGENT!
## number of iterations= 1000
## [1] 44
## [1] "Fitting GMM for chr14 18173523:107043718 iteration 1"
## number of iterations= 525
## [1] "Fitting GMM for chr14 18173523:107043718 iteration 2"
## number of iterations= 547
## [1] "Fitting GMM for chr14 18173523:107043718 iteration 3"
## number of iterations= 392
## [1] "Fitting GMM for chr14 18173523:107043718 iteration 4"
## number of iterations= 429
## [1] 44
## [1] "Fitting GMM for chr15 19725254:101991189 iteration 1"
## WARNING! NOT CONVERGENT!
## number of iterations= 1000
## [1] "Fitting GMM for chr15 19725254:101991189 iteration 2"
## number of iterations= 911
## [1] "Fitting GMM for chr15 19725254:101991189 iteration 3"
## WARNING! NOT CONVERGENT!
## number of iterations= 1000
## [1] "Fitting GMM for chr15 19725254:101991189 iteration 4"
## WARNING! NOT CONVERGENT!
## number of iterations= 1000
## [1] 44
## [1] "Fitting GMM for chr16 0:36311158 iteration 1"
## number of iterations= 862
## [1] "Fitting GMM for chr16 0:36311158 iteration 2"
## WARNING! NOT CONVERGENT!
## number of iterations= 1000
## [1] "Fitting GMM for chr16 0:36311158 iteration 3"
## number of iterations= 747
## [1] "Fitting GMM for chr16 0:36311158 iteration 4"
## number of iterations= 929
## [1] 44
## [1] "Fitting GMM for chr16 36334460:90338345 iteration 1"
## number of iterations= 745
## [1] "Fitting GMM for chr16 36334460:90338345 iteration 2"
## number of iterations= 770
## [1] "Fitting GMM for chr16 36334460:90338345 iteration 3"
## number of iterations= 663
## [1] "Fitting GMM for chr16 36334460:90338345 iteration 4"
## number of iterations= 853
## [1] 44
## [1] "Fitting GMM for chr17 0:22813679 iteration 1"
## WARNING! NOT CONVERGENT!
## number of iterations= 1000
## [1] "Fitting GMM for chr17 0:22813679 iteration 2"
## WARNING! NOT CONVERGENT!
## number of iterations= 1000
## [1] "Fitting GMM for chr17 0:22813679 iteration 3"
## WARNING! NOT CONVERGENT!
## number of iterations= 1000
## [1] "Fitting GMM for chr17 0:22813679 iteration 4"
## WARNING! NOT CONVERGENT!
## number of iterations= 1000
## [1] 44
## [1] "Fitting GMM for chr17 26566633:83257441 iteration 1"
## number of iterations= 900
## [1] "Fitting GMM for chr17 26566633:83257441 iteration 2"
## number of iterations= 652
## [1] "Fitting GMM for chr17 26566633:83257441 iteration 3"
## number of iterations= 838
## [1] "Fitting GMM for chr17 26566633:83257441 iteration 4"
## number of iterations= 661
## [1] 44
## [1] "Fitting GMM for chr18 20861206:80373285 iteration 1"
## WARNING! NOT CONVERGENT!
## number of iterations= 1000
## [1] "Fitting GMM for chr18 20861206:80373285 iteration 2"
## WARNING! NOT CONVERGENT!
## number of iterations= 1000
## [1] "Fitting GMM for chr18 20861206:80373285 iteration 3"
## WARNING! NOT CONVERGENT!
## number of iterations= 1000
## [1] "Fitting GMM for chr18 20861206:80373285 iteration 4"
## number of iterations= 161
## [1] 44
## [1] "Fitting GMM for chr19 0:24498980 iteration 1"
## WARNING! NOT CONVERGENT!
## number of iterations= 1000
## [1] "Fitting GMM for chr19 0:24498980 iteration 2"
## WARNING! NOT CONVERGENT!
## number of iterations= 1000
## [1] "Fitting GMM for chr19 0:24498980 iteration 3"
## number of iterations= 894
## [1] "Fitting GMM for chr19 0:24498980 iteration 4"
## WARNING! NOT CONVERGENT!
## number of iterations= 1000
## [1] 44
## [1] "Fitting GMM for chr19 27190874:58617616 iteration 1"
## WARNING! NOT CONVERGENT!
## number of iterations= 1000
## [1] "Fitting GMM for chr19 27190874:58617616 iteration 2"
## WARNING! NOT CONVERGENT!
## number of iterations= 1000
## [1] "Fitting GMM for chr19 27190874:58617616 iteration 3"
## WARNING! NOT CONVERGENT!
## number of iterations= 1000
## [1] "Fitting GMM for chr19 27190874:58617616 iteration 4"
## WARNING! NOT CONVERGENT!
## number of iterations= 1000
## [1] 44
## [1] "Fitting GMM for chr20 0:26436232 iteration 1"
## number of iterations= 908
## [1] "Fitting GMM for chr20 0:26436232 iteration 2"
## WARNING! NOT CONVERGENT!
## number of iterations= 1000
## [1] "Fitting GMM for chr20 0:26436232 iteration 3"
## WARNING! NOT CONVERGENT!
## number of iterations= 1000
## [1] "Fitting GMM for chr20 0:26436232 iteration 4"
## number of iterations= 846
## [1] 44
## [1] "Fitting GMM for chr20 30038348:64444167 iteration 1"
## number of iterations= 331
## [1] "Fitting GMM for chr20 30038348:64444167 iteration 2"
## number of iterations= 307
## [1] "Fitting GMM for chr20 30038348:64444167 iteration 3"
## number of iterations= 451
## [1] "Fitting GMM for chr20 30038348:64444167 iteration 4"
## number of iterations= 313
## [1] 44
## [1] "Fitting GMM for chr21 12915808:46709983 iteration 1"
## number of iterations= 444
## [1] "Fitting GMM for chr21 12915808:46709983 iteration 2"
## number of iterations= 804
## [1] "Fitting GMM for chr21 12915808:46709983 iteration 3"
## number of iterations= 587
## [1] "Fitting GMM for chr21 12915808:46709983 iteration 4"
## number of iterations= 755
## [1] 44
## [1] "Fitting GMM for chr22 15054318:50818468 iteration 1"
## number of iterations= 657
## [1] "Fitting GMM for chr22 15054318:50818468 iteration 2"
## number of iterations= 786
## [1] "Fitting GMM for chr22 15054318:50818468 iteration 3"
## number of iterations= 720
## [1] "Fitting GMM for chr22 15054318:50818468 iteration 4"
## number of iterations= 971
## [1] 44
After estimating CNVs, we cluster the cells into \(k = 3\) clusters, with the hope of finding one large cluster of normal cells and two smaller clusters composed of CAFs and PDAC cells.
bic_table <- read.table("./Elyada_BIC_LR.txt",
sep = "\t",
row.names = 1,
header = TRUE,
check.names = FALSE)
cand_regions <- rownames(bic_table[bic_table$`BIC difference` > 1000 & bic_table$`LRT adj. p-val` < .01, ])
hist1 <- plotHistogram(cnv_est[, cand_regions],
cpm,
clusters = 3,
zscoreThreshold = 3,
celltypes = pdac$SingleR.labels.bulk,
patients = pdac$sample)We add the normal vs. malignant cell labels in to our Seurat object’s metadata, then visualize the results. As expected, the malignant cells are located in the clusters identified by SingleR as ductal cells and fibroblasts. This indicates that CONICSmat did a solid job of estimating the CNVs - no easy feat with sparse, noisy single-cell data.
normal <- which(hist1 == 2)
malignant <- which(hist1 != 2)
pdac@meta.data$malig <- ifelse(rownames(pdac@meta.data) %in% names(normal), "Normal", "Malignant")
p7 <- DimPlot(pdac, reduction = "tsne", group.by = "malig") +
scale_color_manual(values = wes_palette("Zissou1", n = 5)[c(5, 2)]) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme(plot.title = element_blank()) +
theme_yehlab() +
guides(color = guide_legend(nrow = 1, override.aes = list(size = 3)))
p7## used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
## Ncells 8277599 442.1 51747732 2763.7 NA 31066857 1659.2
## Vcells 269765466 2058.2 1666358438 12713.4 16384 2082948016 15891.7
Next, we use Dr. Xianlu Peng’s DECODER in order to deconvolve the dataset and assign weights to each cell using non-negative matrix factorization. We calculate basal & classical PDAC, normal & activated stroma, immune, and endocrine & exocrine pancreas compartment weights.
sample_wts_unscaled <- Decon_single_sample("TCGA_RNAseq_PAAD", pdac@assays$SCT@data, "geneSymbol")
sample_wts <- Norm_PDAC_weights(sample_wts_unscaled)
pdac <- AddMetaData(pdac, sample_wts$Immune, "immune")
pdac <- AddMetaData(pdac, sample_wts$bcRatio, "bc_ratio")
pdac <- AddMetaData(pdac, sample_wts$Exocrine, "exocrine")
pdac <- AddMetaData(pdac, sample_wts$Endocrine, "endocrine")
pdac <- AddMetaData(pdac, sample_wts_unscaled[, 9], "basal")
pdac <- AddMetaData(pdac, sample_wts_unscaled[, 5], "classical")
pdac <- AddMetaData(pdac, sample_wts$NormalStroma, "norm_stroma")
pdac <- AddMetaData(pdac, sample_wts$ActivatedStroma, "act_stroma")The basal weights are highest in a subcluster of the ductal group identified through SingleR. This is interesting as the authors did not find evidence of the basal subtype in their paper.
p8 <- FeaturePlot(pdac, reduction = "tsne", features = "basal") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme(plot.title = element_blank()) +
theme_yehlab() +
NoLegend()
p8The classical weights are highest in another subcluster of the ductal cluster, and cells with high classical weights do not collocate with those that have high basal weights. The putative classical and basal PDAC cells also align closely with the cells identified through CONCISmat as malignant.
p9 <- FeaturePlot(pdac, reduction = "tsne", features = "classical") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme(plot.title = element_blank()) +
theme_yehlab() +
NoLegend()
p9The cluster identified through SingleR as acinar cells is the only cluster with high exocrine pancreas weights.
p10 <- FeaturePlot(pdac, reduction = "tsne", features = "exocrine") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme(plot.title = element_blank()) +
theme_yehlab() +
NoLegend()
p10No cells have high endocrine pancreas weights.
p11 <- FeaturePlot(pdac, reduction = "tsne", features = "endocrine") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme(plot.title = element_blank()) +
theme_yehlab() +
NoLegend()
p11Once again, we confirm the largeness of the immune cell population in this dataset.
p12 <- FeaturePlot(pdac, reduction = "tsne", features = "immune") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme(plot.title = element_blank()) +
theme_yehlab() +
NoLegend()
p12Cells with high normal stroma stroma weights are located in the cluster identified by SingleR as being composed of endothelial cells.
p13 <- FeaturePlot(pdac, reduction = "tsne", features = "norm_stroma") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme(plot.title = element_blank()) +
theme_yehlab() +
NoLegend()
p13Cells with high activated stroma weights are also located in the fibroblast cluster, and do not intersect with the cells that have high normal stroma scores. This indicates that SCISSORS will most likely perform well on the fibroblast cluster and be able to quickly tease out the cell subtypes.
p14 <- FeaturePlot(pdac, reduction = "tsne", features = "act_stroma") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme(plot.title = element_blank()) +
theme_yehlab() +
NoLegend()
p14Now that we have rough labels from SingleR, CNVs from CONICSmat, and compartment weights from DECODER, we should have more than enough cell-level metadata to look for and annotate cell subtypes using SCISSORS.
The fibroblast marker genes provided by Elyada et al match the SingleR results defining cluster 8 as containing fibroblasts.
p15 <- FeaturePlot(pdac, reduction = "tsne", features = "COL1A1") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(title = "COL1A1") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p16 <- FeaturePlot(pdac, reduction = "tsne", features = "COL3A1") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(title = "COL3A1") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p17 <- FeaturePlot(pdac, reduction = "tsne", features = "LUM") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(title = "LUM") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p18 <- FeaturePlot(pdac, reduction = "tsne", features = "DCN") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(title = "DCN") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
(p15 | p16) / (p17 | p18)Here we use ReclusterCells() to identify subclusters within cluster 7. We find four distinct subclusters.
fibro <- ReclusterCells(pdac,
which.clust = 8,
n.variable.genes = 4000,
n.PC = 20,
resolution.vals = c(.1, .2, .3, .4),
k.vals = c(5, 10, 20, 30),
which.dim.reduc = "tsne",
redo.embedding = TRUE,
do.plot = FALSE)
fibro_pc <- Embeddings(fibro, "pca")## [1] "Reclustering cells in cluster 8 using k = 20 & resolution = 0.1, which achieved silhouette score: 0.536"
We’ll again run Fit-SNE on the reclustered cells, for consistencies sake.
# import data
fibro_pc = r.fibro_pc
# run Fit-SNE
affin_fibro = PerplexityBasedNN(fibro_pc, perplexity=40, metric='cosine', random_state=629)
init = initialization.pca(fibro_pc, random_state=629)
tsne_f = TSNEEmbedding(init, affin_fibro, negative_gradient_method='fft')
embed_f1 = tsne_f.optimize(n_iter=250, exaggeration=5, momentum=0.4)
embed_f2 = embed_f1.optimize(n_iter=750, exaggeration=1, momentum=0.8)Pulling the results back into R and visualizing them shows clear separation between the subclusters. There’s some noise in subcluster 0, but other than that the reembedding looks solid.
embed_fibro <- as.matrix(py$embed_f2)
rownames(embed_fibro) <- colnames(fibro)
fibro@reductions$bh_tsne <- fibro@reductions$tsne
fibro@reductions$tsne<- CreateDimReducObject(embeddings = embed_fibro,
key = "FitSNE_",
assay = "SCT",
global = TRUE)
p19 <- DimPlot(fibro, reduction = "tsne") +
scale_color_manual(values = paletteer_d("miscpalettes::brightPastel")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
guides(color = guide_legend(nrow = 1, override.aes = list(size = 3)))
p19Next, we use the VAM method of single cell gene set enrichment analysis to determine which clusters are enriched for the iCAF and myCAF marker genes, as well as the general pan-CAF marker set. We use the marker genes identified by the authors.
icaf_genes <- c("IL6", "PDGFRA", "CXCL12", "CFD", "LMNA",
"AGTR1", "HAS1", "CXCL1", "CXCL2", "CCL2", "IL8")
mycaf_genes <- c("ACTA2", "TAGLN", "MMP11", "MYL9",
"HOPX", "POSTN", "TPM1", "TPM2")
pan_caf_genes <- c("COL1A1", "FAP", "PDPN", "DCN", "VIM")
gene_sets <- list(icaf_genes, mycaf_genes, pan_caf_genes)
names(gene_sets) <- c("iCAF", "myCAF", "Pan-CAF")
for (i in seq(gene_sets)) {
gene_sets[[i]] <- gene_sets[[i]][gene_sets[[i]] %in% rownames(fibro)]
}
fibro <- vamForSeurat(fibro,
gene.set.collection = gene_sets,
gamma = TRUE)
DefaultAssay(fibro) <- "VAMcdf"We can easily define cluster 0 as the myCAF population, and clusters 1 and 2 as the slightly smaller iCAF population. Interestingly, cluster 3 shows no enrichment whatsoever for either the iCAF or myCAF gene sets.
p20 <- FeaturePlot(fibro, reduction = "tsne", features = "iCAF") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(title = "iCAF Enrichment") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p21 <- FeaturePlot(fibro, reduction = "tsne", features = "myCAF") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(title = "myCAF Enrichment") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
(p20 | p21) / p19Interestingly, we have a small cluster of 22 cells that does not appear to express any of the fibroblast, CAF, perivascular, or endothelial markers. We’ll use differential expression testing to find out what makes it unique, and hopefully annotate it. After performing differential expression analysis, we found that the top 3 markers for cluster 4 are CLU, CD74, and CRYAB. Interestingly, CLU and CD74 were found to be differentially expressed in the apCAF population discovered in the KPC mouse models of CAFs in Elyada et al.
DefaultAssay(fibro) <- "SCT"
fibro_markers <- FindAllMarkers(fibro,
assay = "SCT",
logfc.threshold = 1.5,
test.use = "wilcox",
only.pos = TRUE,
random.seed = 629,
verbose = FALSE)
fibro_markers %>%
filter(cluster == 4) %>%
arrange(desc(avg_logFC))| p_val | avg_logFC | pct.1 | pct.2 | p_val_adj | cluster | gene |
|---|
We re-run GSEA, again using the VAM package and the differentially expressed genes for the apCAF population as defined in Elyada et al (with the mouse gene names converted to HGNC symbols). We can see that the apCAF pathway is strongly enriched in cluster 3 as compared to the other CAF clusters.
apcaf_genes <- c("HLA-DQB1", "CD74", "SAA3P", "SLPI")
gene_sets <- list(icaf_genes, mycaf_genes, apcaf_genes, pan_caf_genes)
names(gene_sets) <- c("iCAF", "myCAF", "apCAF", "Pan-CAF")
for (i in seq(gene_sets)) {
gene_sets[[i]] <- gene_sets[[i]][gene_sets[[i]] %in% rownames(fibro)]
}
fibro <- vamForSeurat(fibro,
gene.set.collection = gene_sets,
gamma = TRUE)
DefaultAssay(fibro) <- "VAMcdf"
p22 <- FeaturePlot(fibro, reduction = "tsne", features = "apCAF") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(title = "apCAF Enrichment") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p22 / p19Finally, we add cell labels to our identified clusters.
fibro$label <- case_when(fibro$seurat_clusters == 0 ~ "myCAF",
fibro$seurat_clusters == 1 ~ "iCAF",
fibro$seurat_clusters == 2 ~ "iCAF",
fibro$seurat_clusters == 3 ~ "apCAF")
p23 <- DimPlot(fibro, reduction = "tsne", group.by = "label") +
scale_color_manual(values = paletteer_d("miscpalettes::brightPastel")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
theme(plot.title = element_blank())
p23We hypothesize that the cluster identified as endothelial cells by SingleR contains our normal stroma population, and potentially some other cell types as it contains visible subclusters.
endo <- ReclusterCells(pdac,
which.clust = 12,
n.variable.genes = 4000,
n.PC = 20,
which.dim.reduc = "tsne",
redo.embedding = TRUE,
random.seed = 629)
endo_pc <- Embeddings(endo, "pca")## [1] "Reclustering cells in cluster 12 using k = 50 & resolution = 0.1, which achieved silhouette score: 0.425"
For consistencies’ sake, we’ll run Fit-SNE on the reclustered object.
# import data
endo_pc = r.endo_pc
# run Fit-SNE
affin_endo = PerplexityBasedNN(endo_pc, perplexity=60, metric='cosine', random_state=629)
init = initialization.pca(endo_pc, random_state=629)
tsne_e = TSNEEmbedding(init, affin_endo, negative_gradient_method='fft')
embed_e1 = tsne_e.optimize(n_iter=250, exaggeration=10, momentum=0.6)
embed_e2 = embed_e1.optimize(n_iter=750, exaggeration=1, momentum=0.8)embed_endo <- as.matrix(py$embed_e2)
rownames(embed_endo) <- colnames(endo)
endo@reductions$bh_tsne <- endo@reductions$tsne
endo@reductions$tsne<- CreateDimReducObject(embeddings = embed_endo,
key = "FitSNE_",
assay = "SCT",
global = TRUE)
p24 <- DimPlot(endo, reduction = "tsne") +
scale_color_manual(values = paletteer_d("miscpalettes::brightPastel")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
guides(color = guide_legend(nrow = 1, override.aes = list(size = 3)))
p24By examining the marker genes provided in the Elyada paper, we define cluster 0 as being composed of perivascular cells.
p25 <- FeaturePlot(endo, reduction = "tsne", features = "IGFBP7") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(title = "IGFBP7") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p26 <- FeaturePlot(endo, reduction = "tsne", features = "ACTA2") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(title = "ACTA2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p27 <- FeaturePlot(endo, reduction = "tsne", features = "RGS5") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(title = "RGS5") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
(p25 | p26 | p27) / p24Next, we use the gene markers PLVAP and VWF (again provided in Elyada et al) to identify cluster 1 as the endothelial cells.
p28 <- FeaturePlot(endo, reduction = "tsne", features = "PLVAP") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(title = "PLVAP") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p29 <- FeaturePlot(endo, reduction = "tsne", features = "VWF") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(title = "VWF") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
(p28 | p29) / p24First we add labels to our two clusters, then we visualize the results.
endo$label <- case_when(endo$seurat_clusters == 0 ~ "Perivascular",
endo$seurat_clusters == 1 ~ "Endothelial")
p30 <- DimPlot(endo, reduction = "tsne", group.by = "label") +
scale_color_manual(values = paletteer_d("miscpalettes::brightPastel")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
theme(plot.title = element_blank())
p30Going back to the main Seurat object, we should have a large population of T and NK cells that warrants further investigation. Using CD3D expression we can easily identify clusters 0, 3, and 7 as the mixed T & NK cells. We already see some good separation, so reclustering the cells should have good results.
p31 <- FeaturePlot(pdac, reduction = "tsne", features = "CD3D") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p31 / p4nkt_tumor <- subset(pdac, subset = seurat_clusters %in% c(0, 3, 7) & condition == "PDAC")Here we run ReclusterCells(), while treating the NK & T cell clusters as one large group. This will hopefully allow use to elucidate T cell subtypes. We use fewer PCs for these cells since the differences between immune cells are subtle, and adding more PCs will most likely only contribute noise to our analyses.
nkt_tumor <- ReclusterCells(nkt_tumor,
which.clust = list(0, 3, 7),
merge.clusters = TRUE,
n.variable.genes = 4000,
n.PC = 15,
k.vals = c(30, 40, 50, 60),
resolution.vals = c(.1, .2, .3),
which.dim.reduc = "tsne",
redo.embedding = TRUE,
random.seed = 629)
nkt_tumor_pc <- Embeddings(nkt_tumor, "pca")## [1] "Reclustering cells in cluster 0 using k = 30 & resolution = 0.3, which achieved silhouette score: 0.455"
Once again we’ll run Fit-SNE on the reclustered cells.
# import data
nkt_tumor_pc = r.nkt_tumor_pc
# run Fit-SNE
affin_nkt_tumor = PerplexityBasedNN(nkt_tumor_pc, perplexity=180, metric='cosine', random_state=629)
init = initialization.pca(nkt_tumor_pc, random_state=629)
tsne_nkt_tumor = TSNEEmbedding(init, affin_nkt_tumor, negative_gradient_method='fft')
embed_nkt_tumor1 = tsne_nkt_tumor.optimize(n_iter=250, exaggeration=4, momentum=0.6)
embed_nkt_tumor2 = embed_nkt_tumor1.optimize(n_iter=750, exaggeration=1, momentum=0.8)
affin_nkt_tumor.set_perplexity(75)
embed_nkt_tumor3 = embed_nkt_tumor2.optimize(n_iter=3, exaggeration=2, momentum=0.1)
embed_nkt_tumor4 = embed_nkt_tumor3.optimize(n_iter=500, exaggeration=1, momentum=0.6)The reembedding isn’t perfect, which we somewhat expected as immune cells are difficult to tell apart based on the transcriptome alone.
embed_nkt_tumor <- as.matrix(py$embed_nkt_tumor4)
rownames(embed_nkt_tumor) <- colnames(nkt_tumor)
nkt_tumor@reductions$bh_tsne <- nkt_tumor@reductions$tsne
nkt_tumor@reductions$tsne<- CreateDimReducObject(embeddings = embed_nkt_tumor,
key = "FitSNE_",
assay = "SCT",
global = TRUE)
p32 <- DimPlot(nkt_tumor, reduction = "tsne") +
scale_color_manual(values = paletteer_d("miscpalettes::brightPastel")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
guides(color = guide_legend(nrow = 1, override.aes = list(size = 3)))
p32CLusters 0 and 1 contain our CD4+ T cells, which we characterize using IL7R as we did in the PBMC3k vignette. It’s somewhat outside of our scope here to determine which subtype each cluster belongs to, so we’ll simply assign both clusters the same label and move on.
p33 <- FeaturePlot(nkt_tumor, reduction = "tsne", features = "IL7R") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p34 <- FeaturePlot(nkt_tumor, reduction = "tsne", features = "CD69") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
(p33 | p34) / p32Cluster 3 contains the regulatory T cells.
p35 <- FeaturePlot(nkt_tumor, reduction = "tsne", features = "IL2RA") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p36 <- FeaturePlot(nkt_tumor, reduction = "tsne", features = "FOXP3") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
(p35 | p36) / p32We can find the proliferating T-regs in cluster 7.
p37 <- FeaturePlot(nkt_tumor, reduction = "tsne", features = "TOP2A") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p37 / p32Mast cells can be identified using TPSAB1 expression in cluster 6.
p38 <- FeaturePlot(nkt_tumor, reduction = "tsne", features = "TPSAB1") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p38 / p32NKG7 and PRF1 show us the NK cells in cluster 5.
p39 <- FeaturePlot(nkt_tumor, reduction = "tsne", features = "NKG7") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p40 <- FeaturePlot(nkt_tumor, reduction = "tsne", features = "PRF1") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
(p39 | p40) / p32We use CD8A and CD2 to reveal the CD8+ T cells within clusters 2 and 4.
p41 <- FeaturePlot(nkt_tumor, reduction = "tsne", features = "CD8A") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p42 <- FeaturePlot(nkt_tumor, reduction = "tsne", features = "CD2") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
(p41 | p42) / p32Lastly, we have cluster 8, which doesn’t highly express our pan-T/NK cell markers CD3D and CD2. It could be a myeloid cell cluster that was mistakenly grouped with the T/NK cells. We’ll start with a Wilcoxon test to determine its differentially expressed genes.
clust8_markers <- FindAllMarkers(nkt_tumor,
test.use = "wilcox",
min.diff.pct = .2,
logfc.threshold = .5,
verbose = FALSE,
random.seed = 629) %>%
filter(cluster == 8, p_val_adj < .05) %>%
arrange(desc(1 - p_val_adj)) %>%
print()## p_val avg_logFC pct.1 pct.2 p_val_adj cluster gene
## SPI1 1.276127e-117 0.6394855 0.500 0.015 1.913808e-113 8 SPI1
## TMEM176B.1 1.951396e-81 0.7020478 0.525 0.025 2.926509e-77 8 TMEM176B
## AIF1 5.975173e-61 1.0496310 0.675 0.060 8.960966e-57 8 AIF1
## MS4A6A 1.077447e-49 0.6118500 0.425 0.028 1.615847e-45 8 MS4A6A
## C1QA 4.494148e-41 1.0151405 0.500 0.048 6.739873e-37 8 C1QA
## CST3.1 4.909324e-41 1.3989825 0.700 0.101 7.362513e-37 8 CST3
## CEBPD 3.107994e-39 0.5168535 0.500 0.048 4.661059e-35 8 CEBPD
## C1QC 4.910297e-37 1.0911696 0.425 0.038 7.363972e-33 8 C1QC
## LST1 5.979516e-36 0.7781232 0.625 0.084 8.967481e-32 8 LST1
## GRN 1.074343e-33 0.5834286 0.600 0.082 1.611192e-29 8 GRN
## HLA-DRB5 9.274134e-33 0.9552109 0.775 0.152 1.390842e-28 8 HLA-DRB5
## HLA-DQA1 9.993726e-32 1.0629547 0.700 0.125 1.498759e-27 8 HLA-DQA1
## HLA-DRB1.2 5.443325e-30 1.5532245 0.975 0.292 8.163354e-26 8 HLA-DRB1
## HLA-DRA.1 1.158221e-29 2.2121810 0.950 0.409 1.736985e-25 8 HLA-DRA
## HLA-DMB 4.386335e-29 0.5408872 0.475 0.058 6.578186e-25 8 HLA-DMB
## TYROBP.2 3.246406e-27 1.1379201 0.700 0.141 4.868635e-23 8 TYROBP
## KLF4 1.527066e-25 0.7371558 0.375 0.042 2.290141e-21 8 KLF4
## C1QB 2.688702e-24 0.9941321 0.400 0.051 4.032246e-20 8 C1QB
## FCER1G.2 1.827197e-23 0.8256689 0.550 0.097 2.740247e-19 8 FCER1G
## HLA-DQA2 2.113476e-23 0.5648281 0.375 0.045 3.169581e-19 8 HLA-DQA2
## HLA-DPB1 6.539129e-23 1.5900322 0.900 0.387 9.806732e-19 8 HLA-DPB1
## CD74.3 7.169626e-23 1.9170476 0.975 0.724 1.075229e-18 8 CD74
## LYZ 2.281872e-22 1.3494472 0.550 0.107 3.422123e-18 8 LYZ
## HLA-DQB1 4.352556e-22 1.1332910 0.700 0.188 6.527528e-18 8 HLA-DQB1
## PLAUR.1 1.321580e-21 0.5098414 0.400 0.056 1.981974e-17 8 PLAUR
## HLA-DMA 1.559863e-20 0.6523222 0.575 0.121 2.339327e-16 8 HLA-DMA
## HLA-DPA1.1 2.353354e-20 1.4657017 0.850 0.368 3.529325e-16 8 HLA-DPA1
## NPC2 7.072958e-18 0.6875523 0.625 0.162 1.060731e-13 8 NPC2
## CD83 2.621481e-16 0.8046446 0.700 0.214 3.931435e-12 8 CD83
## TYMP 1.904906e-12 0.5372419 0.725 0.266 2.856788e-08 8 TYMP
## APOE.1 5.771981e-11 0.6848562 0.300 0.061 8.656241e-07 8 APOE
## CXCL8 6.271403e-11 0.5556395 0.250 0.042 9.405223e-07 8 CXCL8
## LGALS1.2 1.515934e-10 0.7172623 0.725 0.319 2.273447e-06 8 LGALS1
## GLUL.1 2.783779e-10 0.5003458 0.475 0.146 4.174833e-06 8 GLUL
## IER3 1.476494e-09 0.5731977 0.400 0.115 2.214298e-05 8 IER3
## GSTP1 1.868260e-08 0.6851803 0.550 0.240 2.801829e-04 8 GSTP1
## ATF3 2.529892e-08 0.7294797 0.525 0.187 3.794079e-04 8 ATF3
## S100A11.2 2.753010e-08 0.6309440 0.725 0.429 4.128689e-04 8 S100A11
## SAT1 8.290989e-08 0.9010269 0.900 0.662 1.243400e-03 8 SAT1
## PSAP 2.522505e-07 0.5706095 0.625 0.332 3.783000e-03 8 PSAP
## PNRC1.1 3.015767e-06 0.5040474 0.850 0.632 4.522746e-02 8 PNRC1
Interestingly, several intermediate monocyte markers - LYZ, HLA-DRA, CD74, and HLA-DPB1 - are differentially expressed in cluster 8.
We’ll plot some of those markers below.
p43 <- FeaturePlot(nkt_tumor, reduction = "tsne", features = "LYZ") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p44 <- FeaturePlot(nkt_tumor, reduction = "tsne", features = "HLA-DRA") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p45 <- FeaturePlot(nkt_tumor, reduction = "tsne", features = "CD74") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p46 <- FeaturePlot(nkt_tumor, reduction = "tsne", features = "HLA-DPB1") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
((p43 | p44) / (p45 | p46)) / p32We add labels to our T cell Seurat object and visualize the results.
nkt_tumor$label <- case_when(nkt_tumor$seurat_clusters == 0 ~ "CD4+ T",
nkt_tumor$seurat_clusters == 1 ~ "CD4+ T",
nkt_tumor$seurat_clusters == 2 ~ "CD8+ T",
nkt_tumor$seurat_clusters == 3 ~ "T-reg",
nkt_tumor$seurat_clusters == 4 ~ "CD8+ T",
nkt_tumor$seurat_clusters == 5 ~ "NK",
nkt_tumor$seurat_clusters == 6 ~ "Mast",
nkt_tumor$seurat_clusters == 7 ~ "Proliferating T-reg",
nkt_tumor$seurat_clusters == 8 ~ "Intermediate Monocyte")
p47 <- DimPlot(nkt_tumor, reduction = "tsne", group.by = "label") +
scale_color_manual(values = paletteer_d("miscpalettes::brightPastel")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
theme(plot.title = element_blank()) +
guides(color = guide_legend(nrow = 2, override.aes = list(size = 3)))
p47nkt_norm <- subset(pdac, subset = seurat_clusters %in% c(0, 3, 7) & condition == "AdjNorm")nkt_norm <- ReclusterCells(nkt_norm,
which.clust = list(0, 3, 7),
merge.clusters = TRUE,
n.variable.genes = 4000,
n.PC = 10,
k.vals = c(5, 10, 15, 20),
resolution.vals = c(.1, .2, .3),
which.dim.reduc = "tsne",
redo.embedding = TRUE,
random.seed = 629)
nkt_norm_pc <- Embeddings(nkt_norm, "pca")## [1] "Reclustering cells in cluster 0 using k = 10 & resolution = 0.3, which achieved silhouette score: 0.439"
As with the other reclusterings, we’ll run Fit-SNE in order to (hopefully) obtain a better low-dimensional embedding of our cells.
# import data
nkt_norm_pc = r.nkt_norm_pc
# run Fit-SNE
affin_nkt_norm = PerplexityBasedNN(nkt_norm_pc, perplexity=100, metric='cosine', random_state=629)
init = initialization.pca(nkt_norm_pc, random_state=629)
tsne_nkt_norm = TSNEEmbedding(init, affin_nkt_norm, negative_gradient_method='fft')
embed_nkt_norm1 = tsne_nkt_norm.optimize(n_iter=250, exaggeration=12, momentum=0.6)
embed_nkt_norm2 = embed_nkt_norm1.optimize(n_iter=750, exaggeration=1, momentum=0.8)
affin_nkt_norm.set_perplexity(20)
embed_nkt_norm3 = embed_nkt_norm2.optimize(n_iter=20, exaggeration=2, momentum=0.7)
embed_nkt_norm4 = embed_nkt_norm3.optimize(n_iter=500)embed_nkt_norm <- as.matrix(py$embed_nkt_norm4)
rownames(embed_nkt_norm) <- colnames(nkt_norm)
nkt_norm@reductions$bh_tsne <- nkt_norm@reductions$tsne
nkt_norm@reductions$tsne<- CreateDimReducObject(embeddings = embed_nkt_norm,
key = "FitSNE_",
assay = "SCT",
global = TRUE)
p48 <- DimPlot(nkt_norm, reduction = "tsne") +
scale_color_manual(values = paletteer_d("miscpalettes::brightPastel")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
guides(color = guide_legend(nrow = 1, override.aes = list(size = 3)))
p48First we’ll use a Wilcoxon test to determine which genes characterize each cluster.
nkt_norm_markers <- FindAllMarkers(nkt_norm,
logfc.threshold = .5,
min.diff.pct = .2,
verbose = FALSE,
only.pos = TRUE,
random.seed = 629) %>%
filter(p_val_adj < .05)
nkt_norm_markers %>%
group_by(cluster) %>%
top_n(n = 2, wt = avg_logFC)| p_val | avg_logFC | pct.1 | pct.2 | p_val_adj | cluster | gene |
|---|---|---|---|---|---|---|
| 0 | 1.1374536 | 0.805 | 0.283 | 0.0e+00 | 0 | KLRB1 |
| 0 | 0.6289870 | 0.893 | 0.586 | 0.0e+00 | 0 | IL7R |
| 0 | 1.1797617 | 0.936 | 0.309 | 0.0e+00 | 1 | GZMK |
| 0 | 0.5931411 | 0.806 | 0.444 | 0.0e+00 | 1 | CST7 |
| 0 | 0.6538218 | 0.450 | 0.047 | 0.0e+00 | 2 | KLRC1 |
| 0 | 0.9775093 | 0.535 | 0.119 | 0.0e+00 | 2 | GNLY |
| 0 | 1.6999634 | 0.798 | 0.085 | 0.0e+00 | 3 | GZMB |
| 0 | 1.8534674 | 0.685 | 0.144 | 0.0e+00 | 3 | GNLY |
| 0 | 0.7129245 | 0.642 | 0.076 | 0.0e+00 | 4 | SELL |
| 0 | 0.6277972 | 0.593 | 0.102 | 0.0e+00 | 4 | CCR7 |
| 0 | 1.1823501 | 0.663 | 0.041 | 0.0e+00 | 5 | TNFRSF4 |
| 0 | 1.1552110 | 0.894 | 0.470 | 0.0e+00 | 5 | LTB |
| 0 | 1.5419795 | 0.453 | 0.048 | 0.0e+00 | 6 | AREG |
| 0 | 0.8713618 | 0.660 | 0.191 | 0.0e+00 | 6 | TYROBP |
| 0 | 2.0733590 | 1.000 | 0.309 | 0.0e+00 | 7 | TUBA1B |
| 0 | 2.1312287 | 1.000 | 0.423 | 3.4e-06 | 7 | HMGB2 |
We can use IL7R and S100A4 expression to identify the memory CD4+ T cells in cluster 0. IL7R and CCR7 identify the naive CD4+ T cells in the adjacent cluster 4.
p49 <- FeaturePlot(nkt_norm, reduction = "tsne", features = "IL7R") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p50 <- FeaturePlot(nkt_norm, reduction = "tsne", features = "S100A4") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p51 <- FeaturePlot(nkt_norm, reduction = "tsne", features = "CCR7") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
(p50 | p49 | p51) / p48Next we reveal the CD8+ T cells in clusters 1 and 2 with CD8A, as per usual.
p52 <- FeaturePlot(nkt_norm, reduction = "tsne", features = "CD8A") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p52 / p48The T-reg cluster, cluster 5, is identified using TIGIT and FOXP3.
p53 <- FeaturePlot(nkt_norm, reduction = "tsne", features = "TIGIT") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p54 <- FeaturePlot(nkt_norm, reduction = "tsne", features = "FOXP3") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
(p53 | p54) / p48The natural killers can be found in cluster 4 through their expression of PRF1 and NKG7. The NKG7 expression also confirms the identities of the CD8+ T cells we just annotated.
p55 <- FeaturePlot(nkt_norm, reduction = "tsne", features = "PRF1") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p56 <- FeaturePlot(nkt_norm, reduction = "tsne", features = "NKG7") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
(p55 | p56) / p48The tiny proliferating T-reg population in cluster 7 is characterized by TOP2A.
p57 <- FeaturePlot(nkt_norm, reduction = "tsne", features = "TOP2A") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p57 / p48After running another Wilcoxon differential expression test specifically for the last cluster, cluster 6, we can see that several of the myeloid marker genes from Elyada et al are upregulated - CD68, LYZ, CD14, HLA-DRA. LST1 is also significantly differentially expressed; LST1 is involved in monocyte differentiation. Lastly, TYROBP, aka DAP12, is associated with inflammation and is highly expressed in activated monocytes.
FindMarkers(nkt_norm,
ident.1 = 6,
only.pos = TRUE,
min.pct = .3,
verbose = FALSE,
random.seed = 629) %>%
filter(p_val_adj < .05) | p_val | avg_logFC | pct.1 | pct.2 | p_val_adj | |
|---|---|---|---|---|---|
| AREG | 0.0e+00 | 1.5419795 | 0.453 | 0.048 | 0.0000000 |
| CD68 | 0.0e+00 | 0.3393221 | 0.340 | 0.034 | 0.0000000 |
| CCDC88A | 0.0e+00 | 0.2926852 | 0.340 | 0.038 | 0.0000000 |
| GRN | 0.0e+00 | 0.4770964 | 0.453 | 0.075 | 0.0000000 |
| LST1 | 0.0e+00 | 0.7170228 | 0.509 | 0.099 | 0.0000000 |
| TYROBP | 0.0e+00 | 0.8713618 | 0.660 | 0.191 | 0.0000000 |
| LMO4 | 0.0e+00 | 0.3398349 | 0.396 | 0.079 | 0.0000000 |
| LTC4S | 0.0e+00 | 0.4283473 | 0.302 | 0.050 | 0.0000000 |
| FCER1G | 0.0e+00 | 0.5973548 | 0.547 | 0.184 | 0.0000000 |
| IFITM3 | 0.0e+00 | 0.3271327 | 0.321 | 0.066 | 0.0000000 |
| AIF1 | 0.0e+00 | 0.6072226 | 0.453 | 0.137 | 0.0000000 |
| APOC1 | 0.0e+00 | 0.3694114 | 0.340 | 0.079 | 0.0000001 |
| CD14 | 0.0e+00 | 0.4623759 | 0.321 | 0.074 | 0.0000003 |
| FTL | 0.0e+00 | 0.8599523 | 0.981 | 0.953 | 0.0000055 |
| CAPG | 0.0e+00 | 0.4751300 | 0.547 | 0.219 | 0.0000092 |
| HLA-DMA | 0.0e+00 | 0.3317938 | 0.453 | 0.149 | 0.0000126 |
| NPC2 | 0.0e+00 | 0.4810930 | 0.566 | 0.239 | 0.0000236 |
| HLA-DMB | 0.0e+00 | 0.2830020 | 0.340 | 0.093 | 0.0000253 |
| IGHM | 0.0e+00 | 0.5623978 | 0.321 | 0.088 | 0.0000362 |
| CST3 | 0.0e+00 | 0.5899121 | 0.491 | 0.204 | 0.0000621 |
| STMN1 | 0.0e+00 | 0.3356627 | 0.321 | 0.091 | 0.0000767 |
| TKT | 0.0e+00 | 0.2536298 | 0.377 | 0.120 | 0.0001986 |
| ASAH1 | 0.0e+00 | 0.2597505 | 0.377 | 0.123 | 0.0004828 |
| CMTM6 | 2.0e-07 | 0.2544468 | 0.415 | 0.153 | 0.0027496 |
| LGALS1 | 3.0e-07 | 0.5851548 | 0.566 | 0.295 | 0.0033734 |
| HLA-DQB1 | 3.0e-07 | 0.4839694 | 0.491 | 0.238 | 0.0042709 |
| ZFP36L1 | 9.0e-07 | 0.3381036 | 0.604 | 0.312 | 0.0112503 |
| LYZ | 1.4e-06 | 0.7266010 | 0.547 | 0.312 | 0.0172122 |
| SAT1 | 2.2e-06 | 0.5650706 | 0.585 | 0.332 | 0.0272172 |
| HLA-DRA | 2.4e-06 | 0.8682362 | 0.679 | 0.567 | 0.0306470 |
| FTH1 | 2.8e-06 | 0.6358502 | 0.981 | 0.977 | 0.0348959 |
| GLUL | 3.7e-06 | 0.3550762 | 0.340 | 0.131 | 0.0473811 |
We plot some of the intermediate monocyte marker genes. We can say with reasonable confidence that cluster 6 contains intermediate monocytes.
p58 <- FeaturePlot(nkt_norm, reduction = "tsne", features = "LYZ") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p59 <- FeaturePlot(nkt_norm, reduction = "tsne", features = "HLA-DRA") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p60 <- FeaturePlot(nkt_norm, reduction = "tsne", features = "LST1") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p61 <- FeaturePlot(nkt_norm, reduction = "tsne", features = "TYROBP") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
((p58 | p59) / (p60 | p61)) / p48We add final cell labels to our 8 cell clusters.
nkt_norm$label <- case_when(nkt_norm$seurat_clusters == 0 ~ "Memory CD4+ T",
nkt_norm$seurat_clusters == 1 ~ "CD8+ T",
nkt_norm$seurat_clusters == 2 ~ "CD8+ T",
nkt_norm$seurat_clusters == 3 ~ "Naive CD4+ T",
nkt_norm$seurat_clusters == 4 ~ "NK",
nkt_norm$seurat_clusters == 5 ~ "T-reg",
nkt_norm$seurat_clusters == 6 ~ "Intermediate Monocyte",
nkt_norm$seurat_clusters == 7 ~ "Proliferating T-reg")
p62 <- DimPlot(nkt_norm, reduction = "tsne", group.by = "label") +
scale_color_manual(values = paletteer_d("miscpalettes::brightPastel")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
theme(plot.title = element_blank()) +
guides(color = guide_legend(nrow = 2, override.aes = list(size = 3)))
p62We use KRT8 expression to show the ductal cells residing in clusters 5 and 10. We note that some regions of clusters 5 and 10 have no KRT8 expression, which likely means that they are composed of different cell types.
p63 <- FeaturePlot(pdac, reduction = "tsne", features = "KRT8") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p63 / p4We run SCISSORS, then use a Wilcoxon differential expression test to determine potential markers for each cluster.
ductal <- ReclusterCells(pdac,
which.clust = list(5, 10),
merge.clusters = TRUE,
n.variable.genes = 4000,
n.PC = 20,
k.vals = c(10, 20, 30, 40),
resolution.vals = c(.1, .2, .3, .4),
which.dim.reduc = "tsne",
redo.embedding = TRUE,
do.plot = FALSE,
random.seed = 629)
ductal_markers <- FindAllMarkers(ductal,
logfc.threshold = .5,
min.diff.pct = .2,
only.pos = TRUE,
verbose = FALSE,
random.seed = 629) %>%
filter(p_val_adj < .05)
ductal_markers %>%
group_by(cluster) %>%
top_n(n = 3, wt = avg_logFC)## [1] "Reclustering cells in cluster 5 using k = 20 & resolution = 0.1, which achieved silhouette score: 0.532"
| p_val | avg_logFC | pct.1 | pct.2 | p_val_adj | cluster | gene |
|---|---|---|---|---|---|---|
| 0 | 1.6982303 | 0.979 | 0.451 | 0 | 0 | S100P |
| 0 | 1.6562074 | 0.984 | 0.538 | 0 | 0 | IFI27 |
| 0 | 1.3876187 | 0.731 | 0.234 | 0 | 0 | TFF2 |
| 0 | 2.1429484 | 0.908 | 0.166 | 0 | 1 | MMP7 |
| 0 | 2.2650786 | 0.844 | 0.260 | 0 | 1 | CLU |
| 0 | 3.4629760 | 0.866 | 0.389 | 0 | 1 | SPP1 |
| 0 | 1.8261467 | 0.802 | 0.235 | 0 | 2 | HLA-DPB1 |
| 0 | 1.6950322 | 0.740 | 0.204 | 0 | 2 | HLA-DPA1 |
| 0 | 1.8014227 | 0.956 | 0.483 | 0 | 2 | HLA-DRA |
| 0 | 2.8948971 | 0.881 | 0.091 | 0 | 3 | S100A2 |
| 0 | 2.3207514 | 0.435 | 0.010 | 0 | 3 | CST1 |
| 0 | 2.0687281 | 0.701 | 0.112 | 0 | 3 | PLAT |
| 0 | 4.4278494 | 1.000 | 0.073 | 0 | 4 | CELA3A |
| 0 | 4.3544435 | 0.988 | 0.072 | 0 | 4 | CTRB1 |
| 0 | 4.9087421 | 1.000 | 0.116 | 0 | 4 | CTRB2 |
| 0 | 0.8412036 | 0.512 | 0.065 | 0 | 5 | XIST |
| 0 | 1.0563148 | 0.720 | 0.314 | 0 | 5 | GPX2 |
| 0 | 0.7973432 | 0.787 | 0.499 | 0 | 5 | LGALS4 |
| 0 | 2.6843866 | 0.718 | 0.019 | 0 | 6 | CRISP3 |
| 0 | 2.6899477 | 0.846 | 0.103 | 0 | 6 | TCN1 |
| 0 | 3.1343739 | 0.321 | 0.014 | 0 | 6 | SCGB3A1 |
| 0 | 3.5917377 | 0.629 | 0.002 | 0 | 7 | APOA4 |
| 0 | 3.8527263 | 0.671 | 0.007 | 0 | 7 | APOA1 |
| 0 | 2.8721232 | 0.900 | 0.038 | 0 | 7 | ALDOB |
Again, we’ll run Fit-SNE on the reclustered cells.
duct_pc <- Embeddings(ductal, reduction = "pca")# import data
duct_pc = r.duct_pc
# run Fit-SNE
affin_duct = PerplexityBasedNN(duct_pc, perplexity=200, metric='cosine', random_state=629)
init = initialization.pca(duct_pc, random_state=629)
tsne_duct = TSNEEmbedding(init, affin_duct, negative_gradient_method='fft')
embed_d1 = tsne_duct.optimize(n_iter=250, exaggeration=10, momentum=0.6)
embed_d2 = embed_d1.optimize(n_iter=750, exaggeration=1, momentum=0.8)
affin_duct.set_perplexity(30)
embed_d3 = embed_d2.optimize(n_iter=60, exaggeration=2, momentum=0.8)
embed_d4 = embed_d3.optimize(n_iter=500)We pull the results into R, making sure to save the Barnes-Hut t-SNE results in another reduction slot.
embed_duct <- as.matrix(py$embed_d4)
rownames(embed_duct) <- colnames(ductal)
ductal@reductions$bh_tsne <- ductal@reductions$tsne
ductal@reductions$tsne<- CreateDimReducObject(embeddings = embed_duct,
key = "FitSNE_",
assay = "SCT",
global = TRUE)
p64 <- DimPlot(ductal, reduction = "tsne") +
scale_color_manual(values = paletteer_d("miscpalettes::brightPastel")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
guides(color = guide_legend(nrow = 1, override.aes = list(size = 3)))
p64We use ANPEP expression to identify the lipid processing ductal cells in cluster 7.
p65 <- FeaturePlot(ductal, reduction = "tsne", features = "ANPEP") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p65 / p64Expression of SOD3 and CFTR reveals the secretory cells in cluster 1.
p66 <- FeaturePlot(ductal, reduction = "tsne", features = "SOD3") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p67 <- FeaturePlot(ductal, reduction = "tsne", features = "CFTR") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
(p66 | p67) / p64We use TFF1 and TFF2 expression to annotate the classical 1 epithelial cells in clusters 0 and 5. We can also see that the two classical 1 clusters are split by the sample from which the cells originate. Going forward, we’ll simply label both clusters as classical 1.
p68 <- FeaturePlot(ductal, reduction = "tsne", features = "TFF1") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p69 <- FeaturePlot(ductal, reduction = "tsne", features = "TFF2") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p70 <- DimPlot(ductal, reduction = "tsne", group.by = "sample") +
scale_color_manual(values = paletteer_d("miscpalettes::brightPastel")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2", title = "Sample") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
(p68 | p70 | p69) / p64The classical 2 cells are located in cluster 6, as evidenced by their expression of CRISP3.
p71 <- FeaturePlot(ductal, reduction = "tsne", features = "CRISP3") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p71 / p64The basal compartment weights from DECODER are highest in cluster 3, which we will denote as being composed of basal-like PDAC.
p72 <- FeaturePlot(ductal, reduction = "tsne", features = "basal") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2", title = "Basal-like") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p72 / p64As seen earlier in the results of the differential expression text, the top three genes for cluster 2 were HLA-DRA, HLA-DPA1, and HLA-DPB1. This points to cluster 2 being some kind of myeloid cell. We use expression of CD14, CD16 aka FCGR3A, and MS4A7 to identify it as containing CD14+ CD16+ monocytes.
p73 <- FeaturePlot(ductal, reduction = "tsne", features = "CD14") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p74 <- FeaturePlot(ductal, reduction = "tsne", features = "FCGR3A") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p75 <- FeaturePlot(ductal, reduction = "tsne", features = "MS4A7") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
(p73 | p74 | p75) / p64All three genes are significantly differentially expressed in cluster 2.
clust2_markers <- FindMarkers(ductal,
ident.1 = 2,
random.seed = 629,
only.pos = TRUE)
clust2_markers %>%
filter(rownames(clust2_markers) %in% c("CD14", "MS4A7", "FCGR3A")) %>%
arrange(desc(1 - p_val_adj))| p_val | avg_logFC | pct.1 | pct.2 | p_val_adj | |
|---|---|---|---|---|---|
| MS4A7 | 0 | 0.9614771 | 0.485 | 0.014 | 0 |
| CD14 | 0 | 1.2069216 | 0.551 | 0.061 | 0 |
| FCGR3A | 0 | 0.6200975 | 0.300 | 0.023 | 0 |
Lastly, we show that cluster 4 is composed of acinar cells using CTRB2.
p76 <- FeaturePlot(ductal, reduction = "tsne", features = "CTRB2") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p76 / p64We add final cluster labels to our Seurat object and visualize the results.
ductal$label <- case_when(ductal$seurat_clusters == 0 ~ "Classical 1",
ductal$seurat_clusters == 1 ~ "Secretory",
ductal$seurat_clusters == 2 ~ "CD16+ Monocyte",
ductal$seurat_clusters == 3 ~ "Basal-like",
ductal$seurat_clusters == 4 ~ "Acinar",
ductal$seurat_clusters == 5 ~ "Classical 1",
ductal$seurat_clusters == 6 ~ "Classical 2",
ductal$seurat_clusters == 7 ~ "Lipid proc.")
p77 <- DimPlot(ductal, reduction = "tsne", group.by = "label") +
scale_color_manual(values = paletteer_d("miscpalettes::brightPastel")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
theme(plot.title = element_blank()) +
guides(color = guide_legend(nrow = 2, override.aes = list(size = 3)))
p77We’ll run SCISSORS on cluster 11 in order to reveal its subgroups, the plasma cells and plasmacytoid DCs, as defined by their expression of IGJ (aka JCHAIN) and IRF7, respectively.
p78 <- FeaturePlot(pdac, reduction = "tsne", features = "JCHAIN") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2", title = "JCHAIN / IGJ") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p79 <- FeaturePlot(pdac, reduction = "tsne", features = "IRF7") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2", title = "IRF7") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
(p78 | p79) / p4We run SCISSORS with default parameters.
plas <- ReclusterCells(pdac,
which.clust = 11,
n.PC = 20,
which.dim.reduc = "tsne",
redo.embedding = TRUE,
do.plot = FALSE,
random.seed = 629)
plas_pc <- Embeddings(plas, reduction = "pca")## [1] "Reclustering cells in cluster 11 using k = 25 & resolution = 0.1, which achieved silhouette score: 0.632"
Even though in this case it’s probably not necessary, we’ll run Fit-SNE on the cells just so that our axis labels are consistent.
# import data
plas_pc = r.plas_pc
# run Fit-SNE
affin_plas = PerplexityBasedNN(plas_pc, perplexity=30, metric='cosine', random_state=629)
init = initialization.pca(plas_pc, random_state=629)
tsne_plas = TSNEEmbedding(init, affin_plas, negative_gradient_method='fft')
embed_p1 = tsne_plas.optimize(n_iter=250, exaggeration=12, momentum=0.6)
embed_p2 = embed_p1.optimize(n_iter=750, exaggeration=1, momentum=0.8)embed_plas <- as.matrix(py$embed_p2)
rownames(embed_plas) <- colnames(plas)
plas@reductions$bh_tsne <- plas@reductions$tsne
plas@reductions$tsne<- CreateDimReducObject(embeddings = embed_plas,
key = "FitSNE_",
assay = "SCT",
global = TRUE)
p80 <- DimPlot(plas, reduction = "tsne") +
scale_color_manual(values = paletteer_d("miscpalettes::brightPastel")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
guides(color = guide_legend(nrow = 1, override.aes = list(size = 3)))
p80We use the aforementioned markers to identify our two clusters.
p81 <- FeaturePlot(plas, reduction = "tsne", features = "JCHAIN") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2", title = "JCHAIN") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p82 <- FeaturePlot(plas, reduction = "tsne", features = "IRF7") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2", title = "IRF7") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
(p81 | p82) / p80As before, we add cell labels and visualize the results.
plas$label <- case_when(plas$seurat_clusters == 0 ~ "Plasma",
plas$seurat_clusters == 1 ~ "pDC")
p83 <- DimPlot(plas, reduction = "tsne", group.by = "label") +
scale_color_manual(values = paletteer_d("miscpalettes::brightPastel")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
theme(plot.title = element_blank()) +
guides(color = guide_legend(nrow = 1, override.aes = list(size = 3)))
p83We can identify the B cells in cluster 4 using joint expression of MS4A1 and CD79A. The cluster is split into two subclusters by tissue type: adjacent normal and PDAC. While it would be interesting to determine the genetic drivers of that separation, it’s somewhat outside of our scope here.
p84 <- FeaturePlot(pdac, reduction = "tsne", features = "MS4A1") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2", title = "MS4A1") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p85 <- FeaturePlot(pdac, reduction = "tsne", features = "CD79A") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2", title = "CD79A") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p86 <- DimPlot(pdac, reduction = "tsne", group.by = "condition") +
scale_color_manual(values = paletteer_d("miscpalettes::brightPastel")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2", title = "Tissue Type") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
(p84 | p86 | p85) / p4Lastly, we’ll split up the myeloid population by tissue condition (PDAC vs. adjacent normal) just like we did with the NK / T cells. We’ll run SCISSORS, annotate the clusters, and visualize the results.
First we’ll attempt to assign broad cell type labels to each of the four putative myeloid clusters. We’ll use the marker genes from Elyada et al once again.
myo_tumor <- subset(pdac, subset = seurat_clusters %in% c(1, 2, 6, 9) & condition == "PDAC")
p87 <- DimPlot(myo_tumor, reduction = "tsne") +
scale_color_manual(values = paletteer_d("miscpalettes::brightPastel")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
theme(plot.title = element_blank()) +
guides(color = guide_legend(nrow = 1, override.aes = list(size = 3)))
p87Cluster 1 appears to be composed of our resident macrophages due to its expression of CD14 and C1QA.
p88 <- FeaturePlot(myo_tumor, reduction = "tsne", features = "CD14") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2", title = "CD14") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p89 <- FeaturePlot(myo_tumor, reduction = "tsne", features = "C1QA") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2", title = "C1QA") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
(p88 | p89) / p87We can use LYZ and S100A8 expression to reveal the classic monocytes and neutrophils in cluster 2.
p90 <- FeaturePlot(myo_tumor, reduction = "tsne", features = "LYZ") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2", title = "LYZ") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p91 <- FeaturePlot(myo_tumor, reduction = "tsne", features = "S100A8") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2", title = "S100A8") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
(p90 | p91) / p87Cluster 9 seems to house the alternatively activated macrophages, as shown by expression of SPP1.
p92 <- FeaturePlot(myo_tumor, reduction = "tsne", features = "SPP1") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2", title = "SPP1") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p92 / p87Lastly, we show that cluster 6 contains our various DC subtypes through its expression of FCER1A, a canonical dendritic cell marker.
p93 <- FeaturePlot(myo_tumor, reduction = "tsne", features = "FCER1A") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2", title = "FCER1A") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p93 / p87Since we’ve already annotated some of the myeloid cell types, we’ll focus on cluster 2, which contains both monocytes and neutrophils, and the dendritic cells in cluster 6.
myo_reclust <- ReclusterCells(myo_tumor,
which.clust = c(2, 6),
n.PC = 10,
merge.clusters = FALSE,
k.vals = c(10, 20, 30, 40),
resolution.vals = c(.2, .3, .4),
n.variable.genes = 4000,
which.dim.reduc = "tsne",
redo.embedding = TRUE,
do.plot = FALSE,
random.seed = 629)
mono_tumor <- myo_reclust[[1]]
dc_tumor <- myo_reclust[[2]]## [1] "Reclustering cells in cluster 2 using k = 40 & resolution = 0.2, which achieved silhouette score: 0.313"
## [1] "Reclustering cells in cluster 6 using k = 20 & resolution = 0.4, which achieved silhouette score: 0.454"
We’ll again run Fit-SNE on our reclustered cells.
mono_tumor_pc <- Embeddings(mono_tumor, "pca")
dc_tumor_pc <- Embeddings(dc_tumor, "pca")# import data
mono_pc = r.mono_tumor_pc
dc_pc = r.dc_tumor_pc
# Fit-SNE - monocytes
affin_mono = PerplexityBasedNN(mono_pc, perplexity=30, metric='cosine', random_state=629)
init = initialization.pca(mono_pc, random_state=629)
tsne_mono = TSNEEmbedding(init, affin_mono, negative_gradient_method='fft')
embed_mono1 = tsne_mono.optimize(n_iter=250, exaggeration=10, momentum=0.6)
embed_mono2 = embed_mono1.optimize(n_iter=750, exaggeration=1, momentum=0.8)
# Fit-SNE - DC
affin_dc = PerplexityBasedNN(dc_pc, perplexity=60, metric='cosine', random_state=629)
init = initialization.pca(dc_pc, random_state=629)
tsne_dc = TSNEEmbedding(init, affin_dc, negative_gradient_method='fft')
embed_dc1 = tsne_dc.optimize(n_iter=250, exaggeration=8, momentum=0.6)
embed_dc2 = embed_dc1.optimize(n_iter=750, exaggeration=1, momentum=0.8)
affin_dc.set_perplexity(20)
embed_dc3 = embed_dc2.optimize(n_iter=100)We pull the results back into R and visualize them.
embed_mono <- as.matrix(py$embed_mono2)
rownames(embed_mono) <- colnames(mono_tumor)
mono_tumor@reductions$bh_tsne <- mono_tumor@reductions$tsne
mono_tumor@reductions$tsne<- CreateDimReducObject(embeddings = embed_mono,
key = "FitSNE_",
assay = "SCT",
global = TRUE)
p94 <- DimPlot(mono_tumor, reduction = "tsne") +
scale_color_manual(values = paletteer_d("miscpalettes::brightPastel")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
theme(plot.title = element_blank()) +
guides(color = guide_legend(nrow = 1, override.aes = list(size = 3)))
p94embed_dc <- as.matrix(py$embed_dc3)
rownames(embed_dc) <- colnames(dc_tumor)
dc_tumor@reductions$bh_tsne <- dc_tumor@reductions$tsne
dc_tumor@reductions$tsne<- CreateDimReducObject(embeddings = embed_dc,
key = "FitSNE_",
assay = "SCT",
global = TRUE)
p95 <- DimPlot(dc_tumor, reduction = "tsne") +
scale_color_manual(values = paletteer_d("miscpalettes::brightPastel")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
theme(plot.title = element_blank()) +
guides(color = guide_legend(nrow = 1, override.aes = list(size = 3)))
p95First we ID the neutrophils in cluster 1 using S100A8 and S100A9 - marker genes used by Elyada et al.
p96 <- FeaturePlot(mono_tumor, reduction = "tsne", features = "S100A8") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2", title = "S100A8") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p97 <- FeaturePlot(mono_tumor, reduction = "tsne", features = "S100A9") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2", title = "S100A9") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
(p96 | p97) / p94We can identify the classical monocytes in cluster 0 through expression of, you guessed it, CD14, as well as expression of CD16 aka FCGR3A.
p98 <- FeaturePlot(mono_tumor, reduction = "tsne", features = "CD14") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2", title = "CD14") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p99 <- FeaturePlot(mono_tumor, reduction = "tsne", features = "FCGR3A") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2", title = "CD16") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
(p98 | p99) / p94We add cell labels to our Seurat object, then we’re off to the DCs.
mono_tumor$label <- case_when(mono_tumor$seurat_clusters == 0 ~ "Classical Monocyte",
mono_tumor$seurat_clusters == 1 ~ "Neutrophil")We can use CLEC9A to annotate the cDC1 population in cluster 2.
p100 <- FeaturePlot(dc_tumor, reduction = "tsne", features = "CLEC9A") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2", title = "CLEC9A") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p100 / p95Joint expression of CD1A and CD207 reveals the Langerhans-like DCs in cluster 0.
p101 <- FeaturePlot(dc_tumor, reduction = "tsne", features = "CD1A") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2", title = "CD1A") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p102 <- FeaturePlot(dc_tumor, reduction = "tsne", features = "CD207") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2", title = "CD207") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
(p101 | p102) / p95Next we use LAMP3 to identify the activated DCs in cluster 3.
p103 <- FeaturePlot(dc_tumor, reduction = "tsne", features = "LAMP3") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2", title = "LAMP3") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p103 / p95Using a Wilcoxon differential expression test shows that cluster 1 is characterized by several conventionl DC2 marker genes / transcription factors.
dc_tumor_markers <- FindAllMarkers(dc_tumor,
logfc.threshold = .3,
only.pos = TRUE,
verbose = FALSE,
random.seed = 629)
dc_tumor_markers %>%
filter(p_val_adj < .05 & cluster == 1 & gene %in% c("CLEC10A", "KLF4", "ZEB2"))| p_val | avg_logFC | pct.1 | pct.2 | p_val_adj | cluster | gene |
|---|
p104 <- FeaturePlot(dc_tumor, reduction = "tsne", features = "CLEC10A") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p105 <- FeaturePlot(dc_tumor, reduction = "tsne", features = "KLF4") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p106 <- FeaturePlot(dc_tumor, reduction = "tsne", features = "ZEB2") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
(p104 | p105 | p106) / p95We add final subcluster cell type labels to our Seurat object, and then we’re off to the macrophages.
dc_tumor$label <- case_when(dc_tumor$seurat_clusters == 0 ~ "Langerhans-like DC",
dc_tumor$seurat_clusters == 1 ~ "cDC2",
dc_tumor$seurat_clusters == 2 ~ "cDC1",
dc_tumor$seurat_clusters == 3 ~ "Activated DC")p107 <- DimPlot(mono_tumor, reduction = "tsne", group.by = "label") +
scale_color_manual(values = paletteer_d("miscpalettes::brightPastel")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
theme(plot.title = element_blank()) +
guides(color = guide_legend(nrow = 1, override.aes = list(size = 3)))
p107p108 <- DimPlot(dc_tumor, reduction = "tsne", group.by = "label") +
scale_color_manual(values = paletteer_d("miscpalettes::brightPastel")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
theme(plot.title = element_blank()) +
guides(color = guide_legend(nrow = 1, override.aes = list(size = 3)))
p108We have the same four clusters as in the tumor tissue: 1, 2, 6, and 9. We know that clusters 1 and 9 are composed of the resident and alternatively-activated macrophages, respectively. We’ll focus on the monocytes & neutrophils and DCs in clusters 2 and 6, respectively.
myo_norm <- subset(pdac, subset = seurat_clusters %in% c(1, 2, 6, 9) & condition == "AdjNorm")myo_norm_reclust <- ReclusterCells(myo_norm,
which.clust = c(2, 6),
merge.clusters = FALSE,
n.variable.genes = 4000,
n.PC = 10,
which.dim.reduc = "tsne",
redo.embedding = TRUE,
random.seed = 629)
mono_norm <- myo_norm_reclust[[1]]
dc_norm <- myo_norm_reclust[[2]]## [1] "Reclustering cells in cluster 2 using k = 25 & resolution = 0.1, which achieved silhouette score: 0.406"
## [1] "Reclustering cells in cluster 6 using k = 25 & resolution = 0.1, which achieved silhouette score: 0.48"
For the last time, we run Fit-SNE in order to obtain a better embedding.
mono_norm_pc <- Embeddings(mono_norm, "pca")
dc_norm_pc <- Embeddings(dc_norm, "pca")# import data
mono_pc = r.mono_norm_pc
dc_pc = r.dc_norm_pc
# Fit-SNE - monocytes
affin_mono = PerplexityBasedNN(mono_pc, perplexity=30, metric='cosine', random_state=629)
init = initialization.pca(mono_pc, random_state=629)
tsne_mono = TSNEEmbedding(init, affin_mono, negative_gradient_method='fft')
embed_mono1 = tsne_mono.optimize(n_iter=250, exaggeration=10, momentum=0.6)
embed_mono2 = embed_mono1.optimize(n_iter=750, exaggeration=1, momentum=0.8)
# Fit-SNE - DC
affin_dc = PerplexityBasedNN(dc_pc, perplexity=30, metric='cosine', random_state=629)
init = initialization.pca(dc_pc, random_state=629)
tsne_dc = TSNEEmbedding(init, affin_dc, negative_gradient_method='fft')
embed_dc1 = tsne_dc.optimize(n_iter=250, exaggeration=8, momentum=0.6)
embed_dc2 = embed_dc1.optimize(n_iter=750, exaggeration=1, momentum=0.8)We pull the results back into R and visualize them.
embed_mono <- as.matrix(py$embed_mono2)
rownames(embed_mono) <- colnames(mono_norm)
mono_norm@reductions$bh_tsne <- mono_norm@reductions$tsne
mono_norm@reductions$tsne<- CreateDimReducObject(embeddings = embed_mono,
key = "FitSNE_",
assay = "SCT",
global = TRUE)
p109 <- DimPlot(mono_norm, reduction = "tsne") +
scale_color_manual(values = paletteer_d("miscpalettes::brightPastel")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
theme(plot.title = element_blank()) +
guides(color = guide_legend(nrow = 1, override.aes = list(size = 3)))
p109embed_dc <- as.matrix(py$embed_dc2)
rownames(embed_dc) <- colnames(dc_norm)
dc_norm@reductions$bh_tsne <- dc_norm@reductions$tsne
dc_norm@reductions$tsne<- CreateDimReducObject(embeddings = embed_dc,
key = "FitSNE_",
assay = "SCT",
global = TRUE)
p110 <- DimPlot(dc_norm, reduction = "tsne") +
scale_color_manual(values = paletteer_d("miscpalettes::brightPastel")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
theme(plot.title = element_blank()) +
guides(color = guide_legend(nrow = 1, override.aes = list(size = 3)))
p110We use high S100A8 and S100A9 expression to define the neutrophils in cluster 2.
p111 <- FeaturePlot(mono_norm, reduction = "tsne", features = "S100A8") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p112 <- FeaturePlot(mono_norm, reduction = "tsne", features = "S100A9") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
(p111 | p112) / p109Next up are the classical monocytes in clusters 0 and 1, split by sample.
p113 <- FeaturePlot(mono_norm, reduction = "tsne", features = "CD68") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p114 <- FeaturePlot(mono_norm, reduction = "tsne", features = "CD14") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p115 <- DimPlot(mono_norm, reduction = "tsne", group.by = "sample") +
scale_color_manual(values = paletteer_d("miscpalettes::brightPastel")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2", title = "Sample") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
(p113 | p115 | p114) / p109Lastly, it seems we have a small group of CD8+ T cells in cluster 3 that snuck into the myeloid cluster, as defined by their expression of CD2, CD3D (marking them as NK / T cells), and CD8A & NKG7. I don’t believe they’re NK cells as the do not express PRF1 or GZMB.
p116 <- FeaturePlot(mono_norm, reduction = "tsne", features = "CD2") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p117 <- FeaturePlot(mono_norm, reduction = "tsne", features = "CD3D") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p118 <- FeaturePlot(mono_norm, reduction = "tsne", features = "CD8A") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p119 <- FeaturePlot(mono_norm, reduction = "tsne", features = "NKG7") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
((p116 | p117) / (p118 | p119)) / p109We add labels to our clusters.
mono_norm$label <- case_when(mono_norm$seurat_clusters == 0 ~ "Classic Monocyte",
mono_norm$seurat_clusters == 1 ~ "Classic Monocyte",
mono_norm$seurat_clusters == 2 ~ "Neutrophil",
mono_norm$seurat_clusters == 3 ~ "CD8+ T")We annotate cluster 1 as conventional DC1 and cluster 0 as conventional DC2 through mutually exclusive expression of the canonical markers CLEC9A and CLEC10A, respectively.
p120 <- FeaturePlot(dc_norm, reduction = "tsne", features = "CLEC9A") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p121 <- FeaturePlot(dc_norm, reduction = "tsne", features = "CLEC10A") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
(p120 | p121) / p110Labels are added to our adjacent normal tissue DC Seurat object.
dc_norm$label <- case_when(dc_norm$seurat_clusters == 0 ~ "cDC2",
dc_norm$seurat_clusters == 1 ~ "cDC1")Here are the final annotations for the adjacent normal tissue myeloid population.
p122 <- DimPlot(mono_norm, reduction = "tsne", group.by = "label") +
scale_color_manual(values = paletteer_d("awtools::mpalette")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
theme(plot.title = element_blank()) +
guides(color = guide_legend(nrow = 1, override.aes = list(size = 3)))
p122p123 <- DimPlot(dc_norm, reduction = "tsne", group.by = "label") +
scale_color_manual(values = paletteer_d("awtools::mpalette")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
theme(plot.title = element_blank()) +
guides(color = guide_legend(nrow = 1, override.aes = list(size = 3)))
p123We’ll create a quick convenience function to help us save the figures.
saveSCISSORS <- function(plot = NULL,
name = NULL,
border = TRUE,
pub.ready = FALSE) {
if (is.null(plot) | is.null(name)) stop("You forgot some arguments.")
if (pub.ready) {
dir <- "~/Desktop/R/SCISSORS/vignettes/figures_pub/Elyada/"
if (!border) {
plot <- plot +
theme(panel.border = element_blank(),
axis.title = element_blank(),
legend.position = "none")
} else {
plot <- plot +
theme(axis.title = element_blank(),
legend.position = "none")
}
ggsave(filename = paste0(name, ".pdf"),
device = "pdf",
units = "in",
path = dir,
height = 5,
width = 5)
} else {
dir <- "~/Desktop/R/SCISSORS/vignettes/figures_supp/Elyada/"
ggsave(filename = paste0(name, ".pdf"),
device = "pdf",
units = "in",
path = dir,
height = 5,
width = 5)
}
}This section isn’t worth reading; it’s here solely to prove that the figures we present in our publication were dynamically generated during the knitting of this document.
saveSCISSORS(plot = p0, name = "Seurat_Clusters_tSNE.pdf", pub.ready = TRUE, border = FALSE)
saveSCISSORS(plot = p1, name = "Seurat_Clusters_UMAP.pdf", pub.ready = TRUE, border = FALSE)
saveSCISSORS(plot = p4, name = "Seurat_Clusters_FitSNE.pdf", pub.ready = TRUE, border = FALSE)
saveSCISSORS(plot = p5, name = "SingleR_Annos_sc_ref.pdf", pub.ready = TRUE, border = FALSE)
saveSCISSORS(plot = p6, name = "SingleR_Annos_bulk_ref.pdf", pub.ready = TRUE, border = FALSE)
saveSCISSORS(plot = p7, name = "CONICSmat_Annos.pdf", pub.ready = TRUE, border = FALSE)
saveSCISSORS(plot = p8, name = "DECODER_basal.pdf")
saveSCISSORS(plot = p9, name = "DECODER_classical.pdf")
saveSCISSORS(plot = p10, name = "DECODER_exocrine.pdf")
saveSCISSORS(plot = p11, name = "DECODER_endocrine.pdf")
saveSCISSORS(plot = p12, name = "DECODER_immune.pdf")
saveSCISSORS(plot = p13, name = "DECODER_normal_stroma.pdf")
saveSCISSORS(plot = p14, name = "DECODER_act_stroma.pdf")
saveSCISSORS(plot = p15, name = "Fibro_all_cells_COL1A1.pdf")
saveSCISSORS(plot = p16, name = "Fibro_all_cells_COL3A1.pdf")
saveSCISSORS(plot = p17, name = "Fibro_all_cells_LUM.pdf")
saveSCISSORS(plot = p18, name = "Fibro_all_cells_DCN.pdf")
saveSCISSORS(plot = p19, name = "Fibro_SCISSORS_FitSNE.pdf")
saveSCISSORS(plot = p20, name = "Fibro_iCAF_VAM.pdf")
saveSCISSORS(plot = p21, name = "Fibro_myCAF_VAM.pdf")
saveSCISSORS(plot = p22, name = "Fibro_apCAF_VAM.pdf")
saveSCISSORS(plot = p23, name = "Fibro_final_labels.pdf", pub.ready = TRUE, border = FALSE)
saveSCISSORS(plot = p24, name = "Endo_Perivascular_SCISSORS_FitSNE.pdf")
saveSCISSORS(plot = p25, name = "Perivascular_IGFBP7.pdf")
saveSCISSORS(plot = p26, name = "Perivascular_ACTA2.pdf")
saveSCISSORS(plot = p27, name = "Perivascular_RGS5.pdf")
saveSCISSORS(plot = p28, name = "Endothelial_PLVAP.pdf")
saveSCISSORS(plot = p29, name = "Endothelial_VWF.pdf")
saveSCISSORS(plot = p30, name = "Endo_Perivascular_final_labels.pdf", pub.ready = TRUE, border = FALSE)
saveSCISSORS(plot = p31, name = "NKT_all_cells_CD3D.pdf")
saveSCISSORS(plot = p32, name = "NKT_Tumor_SCISSORS_FitSNE.pdf")
saveSCISSORS(plot = p33, name = "NKT_Tumor_CD4T_IL7R.pdf")
saveSCISSORS(plot = p34, name = "NKT_Tumor_CD4T_CD69.pdf")
saveSCISSORS(plot = p35, name = "NKT_Tumor_Treg_IL2RA.pdf")
saveSCISSORS(plot = p36, name = "NKT_Tumor_Treg_FOXP3.pdf")
saveSCISSORS(plot = p37, name = "NKT_Tumor_Prolif_Treg_TOP2A.pdf")
saveSCISSORS(plot = p38, name = "NKT_Tumor_Mast_TPSAB1.pdf")
saveSCISSORS(plot = p39, name = "NKT_Tumor_NK_NKG7.pdf")
saveSCISSORS(plot = p40, name = "NKT_Tumor_NK_PRF1.pdf")
saveSCISSORS(plot = p41, name = "NKT_Tumor_CD8T_CD8A.pdf")
saveSCISSORS(plot = p42, name = "NKT_Tumor_CD8T_CD2.pdf")
saveSCISSORS(plot = p43, name = "NKT_Tumor_Intermediate_Mono_LYZ.pdf")
saveSCISSORS(plot = p44, name = "NKT_Tumor_Intermediate_Mono_HLADRA.pdf")
saveSCISSORS(plot = p45, name = "NKT_Tumor_Intermediate_Mono_CD74.pdf")
saveSCISSORS(plot = p46, name = "NKT_Tumor_Intermediate_Mono_HLADPB1.pdf")
saveSCISSORS(plot = p47, name = "NKT_Tumor_final_labels.pdf", pub.ready = TRUE, border = FALSE)
saveSCISSORS(plot = p48, name = "NKT_Norm_SCISSORS_FitSNE.pdf")
saveSCISSORS(plot = p49, name = "NKT_Norm_CD4T_IL7R.pdf")
saveSCISSORS(plot = p50, name = "NKT_Norm_CD4T_S100A4.pdf")
saveSCISSORS(plot = p51, name = "NKT_Norm_CD4T_CCR7.pdf")
saveSCISSORS(plot = p52, name = "NKT_Norm_CD8T_CD8A.pdf")
saveSCISSORS(plot = p53, name = "NKT_Norm_Treg_TIGIT.pdf")
saveSCISSORS(plot = p54, name = "NKT_Norm_Treg_FOXP3.pdf")
saveSCISSORS(plot = p55, name = "NKT_Norm_NK_PRF1.pdf")
saveSCISSORS(plot = p56, name = "NKT_Norm_NK_NKG7.pdf")
saveSCISSORS(plot = p57, name = "NKT_Norm_Prolif_Treg_TOP2A.pdf")
saveSCISSORS(plot = p58, name = "NKT_Norm_Intermediate_Mono_LYZ.pdf")
saveSCISSORS(plot = p59, name = "NKT_Norm_Intermediate_Mono_HLADRA.pdf")
saveSCISSORS(plot = p60, name = "NKT_Norm_Intermediate_Mono_LST1.pdf")
saveSCISSORS(plot = p61, name = "NKT_Norm_Intermediate_Mono_TYROBP.pdf")
saveSCISSORS(plot = p62, name = "NKT_Norm_SCISSORS_FitSNE.pdf", pub.ready = TRUE, border = FALSE)
saveSCISSORS(plot = p63, name = "Ductal_all_cells_KRT8.pdf")
saveSCISSORS(plot = p64, name = "Ductal_SCISSORS_FitSNE.pdf")
saveSCISSORS(plot = p65, name = "Ductal_Lipid_Proc_ANPEP.pdf")
saveSCISSORS(plot = p66, name = "Ductal_Secretory_SOD3.pdf")
saveSCISSORS(plot = p67, name = "Ductal_Secretory_CFTR.pdf")
saveSCISSORS(plot = p68, name = "Ductal_Classical1_TFF1.pdf")
saveSCISSORS(plot = p69, name = "Ductal_Classical1_TFF2.pdf")
saveSCISSORS(plot = p70, name = "Ductal_Classical1_SampleID.pdf")
saveSCISSORS(plot = p71, name = "Ductal_Classical2_CRISP3.pdf")
saveSCISSORS(plot = p72, name = "Ductal_basal_DECODER.pdf")
saveSCISSORS(plot = p73, name = "Ductal_CD16_Mono_CD14.pdf")
saveSCISSORS(plot = p74, name = "Ductal_CD16_Mono_FCGR3A.pdf")
saveSCISSORS(plot = p75, name = "Ductal_CD16_Mono_MS4A7.pdf")
saveSCISSORS(plot = p76, name = "Ductal_Acinar_CTRB2.pdf")
saveSCISSORS(plot = p77, name = "Ductal_final_labels.pdf", pub.ready = TRUE, border = FALSE)
saveSCISSORS(plot = p78, name = "Plasma_all_cells_IGJ.pdf")
saveSCISSORS(plot = p79, name = "Plasma_all_cells_IRF7.pdf")
saveSCISSORS(plot = p80, name = "Plasma_SCISSORS_FitSNE.pdf")
saveSCISSORS(plot = p81, name = "Plasma_Plasma_JCHAIN.pdf")
saveSCISSORS(plot = p82, name = "Plasma_pDC_IGJ.pdf")
saveSCISSORS(plot = p83, name = "Plasma_final_labels.pdf", pub.ready = TRUE, border = FALSE)
saveSCISSORS(plot = p84, name = "B_all_cells_MS4A1.pdf")
saveSCISSORS(plot = p85, name = "B_all_cells_CD79A.pdf")
saveSCISSORS(plot = p86, name = "B_all_cells_tissue_type.pdf")
saveSCISSORS(plot = p87, name = "Myeloid_Tumor_original_FitSNE.pdf")
saveSCISSORS(plot = p88, name = "Myeloid_Tumor_Resident_Macro_CD14.pdf")
saveSCISSORS(plot = p89, name = "Myeloid_Tumor_Resident_Macro_C1QA.pdf")
saveSCISSORS(plot = p90, name = "Myeloid_Tumor_Classic_Mono_LYZ.pdf")
saveSCISSORS(plot = p91, name = "Myeloid_Tumor_Classic_Mono_S100A8.pdf")
saveSCISSORS(plot = p92, name = "Myeloid_Tumor_Alternative_Macro_Spp1.pdf")
saveSCISSORS(plot = p93, name = "Myeloid_Tumor_DC_FCER1A.pdf")
saveSCISSORS(plot = p94, name = "Monocyte_Tumor_SCISSORS_FitSNE.pdf")
saveSCISSORS(plot = p95, name = "DC_Tumor_SCISSORS_FitSNE.pdf")
saveSCISSORS(plot = p96, name = "Monocyte_Tumor_Neutrophil_S100A8.pdf")
saveSCISSORS(plot = p97, name = "Monocyte_Tumor_Neutrophil_S100A9.pdf")
saveSCISSORS(plot = p98, name = "Monocyte_Tumor_Classic_Mono_CD14.pdf")
saveSCISSORS(plot = p99, name = "Monocyte_Tumor_Classic_Mono_FCGR3A.pdf")
saveSCISSORS(plot = p100, name = "DC_Tumor_cDC1_CLEC9A.pdf")
saveSCISSORS(plot = p101, name = "DC_Tumor_Langerhans_DC_CD1A.pdf")
saveSCISSORS(plot = p102, name = "DC_Tumor_Langerhans_DC_CD207.pdf")
saveSCISSORS(plot = p103, name = "DC_Tumor_Activated_DC_LAMP3.pdf")
saveSCISSORS(plot = p104, name = "DC_Tumor_cDC2_CLEC10A.pdf")
saveSCISSORS(plot = p105, name = "DC_Tumor_cDC2_KLF4.pdf")
saveSCISSORS(plot = p106, name = "DC_Tumor_cDC2_ZEB2.pdf")
saveSCISSORS(plot = p107, name = "Mono_Tumor_final_labels.pdf", pub.ready = TRUE, border = FALSE)
saveSCISSORS(plot = p108, name = "DC_Tumor_final_labels.pdf", pub.ready = TRUE, border = FALSE)
saveSCISSORS(plot = p109, name = "Mono_Norm_SCISSORS_FitSNE.pdf")
saveSCISSORS(plot = p110, name = "DC_Norm_SCISSORS_FitSNE.pdf")
saveSCISSORS(plot = p111, name = "Mono_Norm_Neutrophil_S100A8.pdf")
saveSCISSORS(plot = p112, name = "Mono_Norm_Neutrophil_S100A9.pdf")
saveSCISSORS(plot = p113, name = "Mono_Norm_Classic_Monocyte_CD68.pdf")
saveSCISSORS(plot = p114, name = "Mono_Norm_Classic_Monocyte_CD14.pdf")
saveSCISSORS(plot = p115, name = "Mono_Norm_Classic_MOnocyte_SampleID.pdf")
saveSCISSORS(plot = p116, name = "Mono_Norm_CD8T_CD2.pdf")
saveSCISSORS(plot = p117, name = "Mono_Norm_CD8T_CD3D.pdf")
saveSCISSORS(plot = p118, name = "Mono_Norm_CD8T_CD8A.pdf")
saveSCISSORS(plot = p119, name = "Mono_Norm_CD8T_NKG7.pdf")
saveSCISSORS(plot = p120, name = "DC_Norm_cDC1_CLEC9A.pdf")
saveSCISSORS(plot = p121, name = "DC_Norm_cDC2_CLEC10A.pdf")
saveSCISSORS(plot = p122, name = "Mono_Norm_final_labels.pdf", pub.ready = TRUE, border = FALSE)
saveSCISSORS(plot = p123, name = "DC_Norm_final_labels.pdf", pub.ready = TRUE, border = FALSE)And of course:
sessionInfo()## R version 4.0.3 (2020-10-10)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Catalina 10.15.7
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] parallel stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] wesanderson_0.3.6.9000 reticulate_1.18
## [3] CONICSmat_0.0.0.1 paletteer_1.3.0
## [5] patchwork_1.1.1 mixtools_1.2.0
## [7] SCISSORS_0.0.2.0 SingleCellExperiment_1.12.0
## [9] data.table_1.13.6 cluster_2.1.0
## [11] biomaRt_2.44.1 decoderr_0.0.0.9000
## [13] SingleR_1.2.4 SummarizedExperiment_1.20.0
## [15] Biobase_2.50.0 GenomicRanges_1.42.0
## [17] GenomeInfoDb_1.26.2 IRanges_2.24.1
## [19] S4Vectors_0.28.1 BiocGenerics_0.36.0
## [21] MatrixGenerics_1.2.0 matrixStats_0.57.0
## [23] ggplot2_3.3.3 dplyr_1.0.2
## [25] VAM_0.4.0 MASS_7.3-53
## [27] Seurat_3.2.3 pals_1.6
##
## loaded via a namespace (and not attached):
## [1] tidyselect_1.1.0 RSQLite_2.2.2
## [3] AnnotationDbi_1.50.3 htmlwidgets_1.5.3
## [5] grid_4.0.3 BiocParallel_1.22.0
## [7] Rtsne_0.15 munsell_0.5.0
## [9] codetools_0.2-18 ica_1.0-2
## [11] statmod_1.4.35 scran_1.16.0
## [13] future_1.21.0 miniUI_0.1.1.1
## [15] withr_2.3.0 colorspace_2.0-0
## [17] highr_0.8 knitr_1.30
## [19] ROCR_1.0-11 tensor_1.5
## [21] listenv_0.8.0 labeling_0.4.2
## [23] GenomeInfoDbData_1.2.4 polyclip_1.10-0
## [25] pheatmap_1.0.12 bit64_4.0.5
## [27] farver_2.0.3 parallelly_1.23.0
## [29] vctrs_0.3.6 generics_0.1.0
## [31] xfun_0.20 BiocFileCache_1.12.1
## [33] squash_1.0.9 R6_2.5.0
## [35] ggbeeswarm_0.6.0 rsvd_1.0.3
## [37] locfit_1.5-9.4 bitops_1.0-6
## [39] spatstat.utils_1.20-2 DelayedArray_0.16.0
## [41] assertthat_0.2.1 promises_1.1.1
## [43] scales_1.1.1 beeswarm_0.2.3
## [45] gtable_0.3.0 globals_0.14.0
## [47] goftest_1.2-2 rlang_0.4.10
## [49] splines_4.0.3 lazyeval_0.2.2
## [51] dichromat_2.0-0 prismatic_1.0.0
## [53] BiocManager_1.30.10 yaml_2.2.1
## [55] reshape2_1.4.4 abind_1.4-5
## [57] httpuv_1.5.4 tools_4.0.3
## [59] ellipsis_0.3.1 RColorBrewer_1.1-2
## [61] ggridges_0.5.3 Rcpp_1.0.5
## [63] plyr_1.8.6 progress_1.2.2
## [65] zlibbioc_1.36.0 purrr_0.3.4
## [67] RCurl_1.98-1.2 prettyunits_1.1.1
## [69] rpart_4.1-15 openssl_1.4.3
## [71] deldir_0.2-3 viridis_0.5.1
## [73] pbapply_1.4-3 cowplot_1.1.1
## [75] zoo_1.8-8 ggrepel_0.9.0
## [77] magrittr_2.0.1 RSpectra_0.16-0
## [79] scattermore_0.7 lmtest_0.9-38
## [81] RANN_2.6.1 fitdistrplus_1.1-3
## [83] hms_0.5.3 mime_0.9
## [85] evaluate_0.14 xtable_1.8-4
## [87] XML_3.99-0.5 gridExtra_2.3
## [89] scater_1.16.2 compiler_4.0.3
## [91] tibble_3.0.4 maps_3.3.0
## [93] KernSmooth_2.23-18 crayon_1.3.4
## [95] htmltools_0.5.0 segmented_1.3-1
## [97] mgcv_1.8-33 later_1.1.0.1
## [99] tidyr_1.1.2 DBI_1.1.0
## [101] ExperimentHub_1.14.2 dbplyr_2.0.0
## [103] rappdirs_0.3.1 Matrix_1.3-2
## [105] igraph_1.2.6 pkgconfig_2.0.3
## [107] plotly_4.9.3 vipor_0.4.5
## [109] dqrng_0.2.1 XVector_0.30.0
## [111] stringr_1.4.0 digest_0.6.27
## [113] sctransform_0.3.2 RcppAnnoy_0.0.18
## [115] spatstat.data_1.7-0 rmarkdown_2.6
## [117] leiden_0.3.6 uwot_0.1.10
## [119] edgeR_3.30.3 DelayedMatrixStats_1.10.1
## [121] curl_4.3 kernlab_0.9-29
## [123] shiny_1.5.0 lifecycle_0.2.0
## [125] nlme_3.1-151 jsonlite_1.7.2
## [127] BiocNeighbors_1.6.0 mapproj_1.2.7
## [129] viridisLite_0.3.0 askpass_1.1
## [131] limma_3.44.3 pillar_1.4.7
## [133] lattice_0.20-41 fastmap_1.0.1
## [135] httr_1.4.2 survival_3.2-7
## [137] interactiveDisplayBase_1.26.3 glue_1.4.2
## [139] spatstat_1.64-1 png_0.1-7
## [141] BiocVersion_3.11.1 bit_4.0.4
## [143] nnls_1.4 stringi_1.5.3
## [145] rematch2_2.1.2 blob_1.2.1
## [147] BiocSingular_1.4.0 AnnotationHub_2.20.2
## [149] memoise_1.1.0 irlba_2.3.3
## [151] future.apply_1.7.0